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The  exponential  increase  in  computing  power  that  began  in  the  ’60s  has  allowed  new 
methods  to  evolve  that  would  make  time  dependent  dynamical  studies  possible.  The 
theoretical  calculation  of  cross  sections  of  chemical  reaction  provides  an  understanding 
of  the  details  at  the  molecular  level.  While  experimental  studies  yield  important  data, 
a significant  fraction  of  the  interesting  elementary  reaction  in  nature  involve  short-lived 
molecular  species  that  are  difficult  to  isolate  and  study  in  the  laboratory.  The  prediction 
of  the  rates  and  cross  sections  of  such  chemical  reactions  is  a primary  goal  of  computa- 
tional chemistry  and  molecular  physics.  These  data,  in  addition  to  its  scientific  interest, 
are  important  for  the  modeling  of  complex  chemical  systems  such  as  combustion,  plas- 
mas, and  the  atmosphere. 

In  this  dissertation  we  present  a method  developed  to  take  advantage  of  current  theo- 
retical techniques  by  calculating  quantum  mechanically  resolved  results  using  classical 
nuclei.  Utilizing  the  nearly  classical  behavior  of  nuclei  in  chemical  processes  we  were 
able  to  devise  a system  to  map  quasi-classical  coherent  states  to  trajectories  calculated 


v 


from  the  semi-classical  approximation  of  the  Electron  Nuclear  Dynamics  theory.  Us- 
ing this  new  method  allows  us  to  analyze  the  results  of  Electron  Nuclear  Dynamics 
molecular  scattering  trajectories,  calculated  using  our  computer  code  called  ENDyne, 
and  extract  state  resolved  results.  In  this  way  we  are  able  to  compare  our  results  to 
experimental  data  and  reveal  their  underlying  chemical  and  mechanical  processes.  We 
present  just  such  applications  to  the  vibrationally  resolved  differential  cross  sections  of 
the  reaction  H+  + H20  and  were  able  to  determine  the  mechanics  of  the  non-adiabatic 
inelastic  scattering  channel.  Additionally,  this  result  served  to  verify  the  ability  of  the 
Electron  Nuclear  Dynamics  theory  to  investigate  non-adiabatic  molecular  scattering. 


V 
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CHAPTER  1 
INTRODUCTION 


Physics  cannot  be  reduced  to  a few  trite  formulae.  The  complexity  of  real  world  pro- 
cesses quickly  overwhelm  the  security  we  find  in  fundamental  equations.  Schrodinger’s 
equation  is  one  possible  beginning  of  understanding  molecular  systems  and  has  not  put  a 
single  quantum  chemist  on  unemployment.  In  fact,  early  quantum  mechanics  struggled 
with  applying  this  equation  to  even  the  simplest  chemical  processes.  Without  the  benefit 
of  computers,  novel  solutions  and  approximations  arose  to  apply  a new  understanding 
of  nature  to  chemistry  and  molecular  scattering  dynamics.  Time-dependent  problems 
were  difficult  to  solve  and  even  then  they  were  difficult  to  calculate  quantitatively.  A 
whole  industry  of  activity  arose  around  solving  time-independent  problems  and  devis- 
ing schemes  to  solve  dynamical  problems  to  some  degree  or  another.  Many  of  these 
methods  suffered  from  an  incomplete  description  of  the  dynamical  processes  and  were 
chiefly  restricted  to  adiabatic  processes. 

The  exponential  increase  in  computing  power  that  began  in  the  ’60s  has  allowed 
new  methods  to  evolve  that  make  time-dependent  dynamical  studies  tractable.  The  ab 
initio  calculation  of  cross  sections  of  chemical  reactions  provides  an  understanding  of 
the  details  at  the  molecular  level.  While  experimental  studies  yield  important  data,  a 
significant  fraction  of  the  interesting  elementary  reactions  in  nature  involve  short-lived 
systems  that  are  difficult  to  isolate  and  study  in  the  laboratory.  Prediction  of  the  rates 
and  cross  sections  of  such  chemical  reactions  is  one  primary  goal  of  computational 
chemistry  and  molecular  physics.  These  data,  in  addition  to  its  scientific  interest,  are 
important  for  modeling  complex  chemical  systems  such  as  combustion,  plasmas,  and 
the  atmosphere. 
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Though  highly  desirable,  theoretical  analysis  of  even  the  simplest  cases  of  these 
dynamical  experiments  can  pose  serious  roadblocks  to  the  scientist.  The  numerical  so- 
lutions can  use  massive  amounts  of  computer  time  and  produce  overwhelming  amounts 
of  data.  Theoretical  molecular  dynamics  strives  to  identify  important  parameters  and 
mechanisms  involved  in  the  chemical  or  collision  processes  for  analysis  and  prediction. 
With  this  in  mind,  a process  of  identifying  and  refining  theoretical  methods  leads  to  the 
creation  of  efficient  tools  for  the  study  of  chemical  systems. 

In  our  group  at  the  University  of  Florida,  lead  by  Prof.  N.  Yngve  Ohm  and  Dr.  E. 
Deumens,  we  have  striven  to  create  new  and  efficient  methods  for  both  analysis  and 
prediction  to  apply  to  difficult  problems  in  reactive  dynamics.  Electron  Nuclear  Dy- 
namics (END)  theory  developed  here  provides  approximate  yet  accurate  solutions  to 
time-dependent  problems.  Using  the  time-dependent  variational  principle  (TDVP)  al- 
lows us  to  simultaneously  calculate  both  the  nuclear  and  the  electronic  dynamics  with- 
out the  use  of  predetermined  potential  energy  surfaces  (PES).  Currently  this  theory  is 
implemented  in  the  ENDyne  code  at  its  minimal  level  of  approximation  described  as 
quasi-classical  single-determinantal  electron  dynamics  (QCSD  END).  At  this  level,  us- 
ing classical  nuclei  and  electrons  described  by  a single-determinantal  wave  function,  we 
have  had  success  in  studying  nonadiabatic  dynamics. 

While  we  continue  to  develop  the  END  theory  and  improve  its  implementation  we 
must  pause  and  consider  whether  a more  complex  theory  or  using  more  computer  time 
is  the  only  correct  way  of  pushing  forward.  This  reflection  brings  us  to  the  overrid- 
ing principle  of  this  dissertation:  Using  coherent  states  and  a posteriori  analysis  we 
can  extract  state-resolved  results  from  a semiclassical  description  using  classical  nuclei. 
Specifically,  we  strove  to  study  rovibrational  dynamics  and  the  state-resolved  differ- 
ential cross  sections  in  molecular  scattering.  Building  a technique  on  three  separate 
pillars — Coherent  States,  analysis  of  classical  dynamics,  and  molecular  trajectories  cal- 
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culated  using  QCSD  END  theory — we  successfully  studied  state  resolved  dynamics  of 
several  interesting  systems.  Here  we  present  this  method  and  its  results. 

This  dissertation  is  organized  in  the  following  way.  Chapter  2 introduces  the  scat- 
tering process  we  study.  Time-independent  and  time-dependent  quantum  methods  are 
discussed.  Traditional  semiclassical  methods  are  explained  and  their  relationship  to 
END  is  elucidated.  Chapter  3 is  a review  of  the  END  theory  and  its  implementation. 
Chapter  4,  rovibrational  dynamics,  both  classical  and  quantum,  are  reviewed  in  the  light 
of  scattering  theory.  We  turn  our  attention  to  coherent  states  for  vibrations  and  rota- 
tion in  chapter  5.  The  methods  developed  to  study  vibrational  analysis  are  presented  in 
Chapter  7 and  its  application  to  specific  systems  are  the  subject  of  chapter  8.  Chapters 
9 and  10  does  the  same  for  rotational  analysis.  Finally  we  will  have  some  concluding 
remarks  in  Chapter  1 1 . 


CHAPTER  2 

OVERVIEW  OF  SCATTERING  THEORY 


Scattering  experiments  and  theory  have  been  a cornerstone  of  the  development  of  mod- 
em physics.  From  the  experiments  of  Rutherford  to  modem  high-energy  particle  ac- 
celerators, the  use  of  particle  collisions  has  allowed  scientists  to  infer  the  structure  and 
properties  of  solids,  gases,  molecules,  atoms  and  subatomic  particles.  Well  designed 
and  extremely  precise  experiments  have  lead  to  a wealth  of  data.  Unfortunately,  most 
scattering  experiments  involve  many  collisions  and  complex  dynamics  leading  to  ex- 
perimental observable  that  are  stochastic  in  nature.  Furthermore,  any  experimental  data 
obtained  are  limited  to  an  overall  change  in  the  state  of  the  system  before  and  after  the 
event.  For  the  theorists  this  poses  a two-fold  challenge  to  explain  the  physics  of  the  indi- 
vidual collision  processes  and  then  to  corroborate  those  findings  with  the  experimental 
data. 

The  theoretical  quantity  for  comparison  with  experiment  is  the  cross  section,  al3  ( E ), 
the  notional  area  through  which  a representative  incident  particle  must  pass  at  a given  en- 
ergy to  achieve  a given  change  * — » j in  the  system.  In  molecular  systems  and  processes 
of  chemical  interest  one  must  recognize  that  the  breakdown  of  classical  mechanics  in 
these  regimes  demands  the  use  of  quantum  mechanics.  Quantum  mechanical  effects 
may  come  from  interference,  effects  of  the  uncertainty  principle  and  the  quantum  nature 
of  many  chemical  processes.  On  the  molecular  scale,  an  appropriate  theoretical  basis  is 
that  of  a semiclassical  theory  allowing  the  use  of  our  experience  with  classical  mechan- 
ics and  including  the  relevant  quantum  mechanical  effects.  Our  goal  in  this  chapter  is 
to  present  the  essence  of  molecular  collision  theory  and  its  use  as  an  introduction  to  our 
own  theoretical  developments  in  studying  state  resolved  reaction  dynamics  and  scatter- 
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ing.  After  the  discussion  of  Child[l]  and  Sakurai[2]  we  separate  the  theory  into  three 
broad  classes.  The  three  are  termed: 

Elastic.  Used  in  cases  of  simple  momentum  transfer  for  studying  mean  intermolec- 
ular  potential  functions  averaged  over  internal  states. 

Inelastic.  Used  if  there  is  also  a change  in  internal  energy  yielding  dependence  of 
the  potential  functions  on  the  internal  coordinates  corresponding  to  that  state. 

Reactive.  Used  when  scattering  is  accompanied  by  the  change  in  chemical  compo- 
sition which  reveals  the  nature  of  the  potential  energy  surface  in  the  neighborhood  of 
the  reaction  coordinate. 


The  simplest  molecular  collision  experiments  involve  two  beams  with  defined  veloc- 
ities vi  and  v2  and  a detector  configured  to  observe  a given  transition  from  state  i — >■  j 
of  the  system.  Considerable  simplification  of  subsequent  calculations  can  be  made  by 
transforming  from  the  laboratory  frame  to  the  center  of  mass  coordinate  system  (CM). 

Transforming  the  positions  f\,  r2  and  velocities  v\,  v2  of  the  collision  particles,  with 

— ^ — * 

mass  mi  and  m2  respectively,  to  the  CM  coordinate,  R,  velocity,  V and  internal  coordi- 
nate, r,  velocity,  v; 


2.1  Theoretical  Considerations 


R 


m2  „ 


r 


n - r2 


(2.1) 


and 


v 


Vi  - V2. 


(2.2) 


6 


Then  the  classical  kinetic  energy  T and  the  angular  momentum  L may  be  written 

1111 

T = -mivj  + -m2v%  = -MV2  + -mv2 , 

L = m\  (f*i  x iq)  + m2  (r2  x v2)  = M (^R  x V^j  + m (r  x v) , (2.3) 

where 

„ mim2  .. 

M = mi+m2,  m = . (2.4) 

m\  + m2 

With  the  potential  energy  independent  of  R,  and  of  the  absolute  orientation  of  the 
system,  the  center  of  mass  motion  may  be  factored  out  of  all  equations.  The  remaining 
problem  that  we  confine  ourselves  to  is  equivalent  to  a single  particle  with  reduced  mass 
777,  in  a given  initial  internal  state  i moving  with  translational  energy,  \mv2 , and  orbital 
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angular  momentum  L about  a scattering  center  at  r = 0.  The  inverse  transformation 
between  the  laboratory  and  center  of  mass  frames  is  easily  derived  from  the  preceding 
equations. 

The  collisions  of  molecules  leads  to  scattering  in  all  directions  and  the  resulting 
intensity  distribution  in  the  CM  frame  gives  the  differential  cross  section  (DCS)  defined 
by  the  ratio 


da  ij 


Number  of  scattered  particles  per  unit  time  per  unit  solid  angle 
Number  of  incident  particles  per  unit  time  per  unit  area 


(2.5) 


The  DCS  is  the  fundamental  observable  quantity  in  scattering  experiments  to  which  all 
other  quantities  may  be  related.  The  first  derived  relevant  quantity  is  the  integral  cross 
section 

r 2n  P7T 

<Jij  ( E ) = / / Iij  (9,  </>)  sin  9 d9d(J).  (2.6) 

J o J 0 

In  the  case  when  the  investigation  of  a bulk  phase  is  desired  the  rate  constant  can  be 
related  to  the  cross  section  by 


k (T)  = (kij  (T)}  = v<Tij  P {v)  v 2 dv) 


(2.7) 


where  P ( v ) denotes  the  relative  velocity  distribution,  and  the  brackets  indicate  an  av- 
erage over  initial  internal  states.  In  the  case  of  a gas  transport  processes  the  thermal 
conductivity  depends  on  the  collision  integral 


(2kT\~2  r 

a,  = — / 

\ 7T  m J J o 


/ mv 
\2kT 


(v) exp  - 


2 kTJ 


dv, 


(2.8) 
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where  ov  (v)  is  the  weighted  cross  section 


Jo  Jo 


(2.9) 


In  such  bulk  phase  measurements  much  of  the  details  of  individual  collision  processes 
is  lost  in  the  averaging  involved  in  the  observations.  The  goal  of  molecular  collision 
theory  according  to  Child[l]  is  to  identify  the  dominant  characteristics  of  a given  model 
potential  function  and  to  extract  from  them  all  possible  observable  consequences. 


Classical  mechanics  provides  the  roots  for  modem  physics.  It  provides  a familiar 
basis  to  begin  a more  intricate  investigation  but  also  quantum  mechanics  still  depends 
on  classical  physics  through  the  construction  of  Hamiltonians  and  Lagrangians.  The 
purpose  of  this  section  is  to  provide  a foundation  for  our  investigation  into  QM  scatter- 
ing. The  discussion  leads  from  classical  trajectory  to  the  collision  cross  sections,  and 
concludes  with  some  comments  about  their  validity. 

In  the  CM  system  the  conservation  of  energy  E and  angular  momentum  L given  by 

Eqn.  2.3  gives  the  conditions  of  a given  event  defined  by  the  incident  relative  velocity 

— * 

v and  the  impact  parameter  b in  Fig.  2.2.  The  trajectory  must  lie  in  a plane  since  L is 
conserved.  The  equations  of  motion  may  be  written  in  polar  coordinates: 


2.2  Classical  Theory  of  Scattering 


(2.10) 


(2.11) 


9 


Elimination  of  the  time  derivatives  gives 

^=(^=± (2.12) 

dr  \r ) r2[l-b2/r2  -V(r)/E]2 

the  sign  given  by  the  radial  velocity;  positive  for  outward  and  negative  for  inward  mo- 
tion, respectively. 


The  complete  classical  trajectory  may  now  be  determined  by  integrating  Eqn.  2.12 
from  oo  to  a and  out  again.  For  our  interests,  only  the  overall  effects  of  the  collision  are 
observed  after  the  event.  The  relevant  quantity  is  termed  the  deflection  function  (DF), 
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0,  related  to  E and  L by 


&{E,L)=  7T- 


(2.13) 


which  gives  the  functional  dependence  to  calculate  the  deflection  angle  of  a single  col- 
lision trajectory. 

We  now  move  our  focus  from  a single  trajectory  specified  by  an  initial  energy  and 
impact  parameter  to  reflect  the  inability  of  experiments  on  a molecular  scale  to  specify 
a single  impact  parameter.  As  we  have  discussed,  the  experimental  results  are  reported 
in  terms  of  the  differential  cross  section,  dcr^/dQ,  and  the  collision  cross  section,  oy, 
which  are  notional  areas  leading  to  a particular  angle  and/or  state.  In  experiments  it 
is  not  possible  to  distinguish  between  positive  and  negative  deflection  angles  so  we 
introduce  the  scattering  angle,  9 = |0|,  which  is  the  magnitude  of  the  DF  and  is  limited 
by  0 < 9 < 7r.  We  observe  that  apart  from  discontinuities  at  9 = 0,  9 varies  smoothly 
with  b,  and  hence  trajectories  within  an  impact  parameter  range  from  b to  b + db  defining 
an  area  2nb  d b,  are  deflected  into  an  appropriate  solid  angle  df)  = 2tt  sin  9 d 9.  The  DCS 
is  therefore 


If  more  than  one  value  of  b leads  to  the  same  angle  we  can  appeal  to  the  experiment 
and  observe  that  the  classical  intensity  I is  increased.  Then  it  can  be  replaced  by  an 
appropriate  sum 


(2.14) 


dfl  Y sin9d9/dbk 


(2.15) 
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where  k is  an  index  for  all  values  of  b leading  to  the  same  deflection  angle  9.  The  total 
elastic  cross  section  is  then  defined  as  t he  integral  of  the  differential  form: 

a(E)=  /^-dfi  = 2tt  [*  1(6,  E)  sin6 dO.  (2.16) 

J dlz  Jo 

We  observe  that  the  expression  for  ^ displays  two  possible  singularities,  the  so- 
called  glory  and  rainbow  effects  in  the  classical  differential  cross  section.  The  glory 
effect,  which  may  occur  at  9 = 0 (forward)  or  9 = 7r  (backward),  is  due  to  the  sin  9 in 
the  denominator  evaluating  to  zero  in  Eqn.  2.15.  The  rainbow  singularity,  on  the  other 
hand  corresponds  to  a extremum  in  the  deflection  function,  at  which  d9/db  = 0.  At  this 
point  the  density  of  classical  trajectories  rises  to  infinity  in  analogy  to  the  phenomenon 
in  optics  which  leads  to  rainbows. 

The  principle  deficiency  of  the  classical  theory  is  of  course  that  no  particle  trajectory 
can  ever  be  precisely  defined,  because  of  the  inherent  uncertainties  associated  with  it. 
The  fascinating  effect  of  interference  reveals  the  failure  of  classical  mechanics  at  the 
rainbow  angle  where  many  trajectories  contribute  to  eliminate  the  singularity. 

2.3  Quantum  Mechanical  Scattering  Theory 

The  standard  treatment  of  quantum  mechanical  scattering  theory  employs  a station- 
ary wave  function  and  its  associated  interference  pattern  in  place  of  the  limiting  deflec- 
tion for  a classical  trajectory.  The  wave  functions  represents  the  whole  system  and  is 
independent  of  time  describing  the  behavior  of  the  system  over  a long  period.  This  is  in 
contrast  to  time  dependent  scattering  theory  which  we  will  discuss  later  in  this  chapter. 
The  starting  point  for  this  formulation  is  determining  the  asymptotic  form  of  the  time- 
independent  scattered  wave  function  and  relate  its  angular  dependence  to  the  effect  of 
the  scattering  potential. 
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Formal  scattering  theory  has  been  in  development  since  the  beginning  of  the  formu- 
lation of  quantum  mechanics  until  the  1970’s[3,  4].  The  main  concern  of  the  discipline 
is  to  rigorously  define  a scattering  process  in  mathematical  terms,  to  attempt  a formal, 
exact  solution  of  the  time-dependent  and  the  time-independent  Schrodinger  equations, 
and  to  define  the  formal  tools  to  treat  a scattering  process.  The  purely  formal  solutions 
are  usually  obtained  from  an  integral  equation  reformulation  of  the  original  Schrodinger 
equation  in  both  the  time-dependent  and  the  time-independent  cases.  These  develop- 
ments are  of  high  theoretical  value  but  of  little  practical  consequences  since  no  definite 
algorithm  to  solve  the  scattering  problem  can  be  devised  in  this  way.  Therefore,  a com- 
plete presentation  of  the  formal  theory  is  beyond  the  scope  of  our  discussion.  Rather, 
the  main  concepts  of  the  formal  theory  are  discussed  along  with  an  introduction  to  the 
practical  methods. 

2.3.1  Time-Independent  vs.  Time-Dependent  Scattering 

A scattering  process  is  a real  time-dependent  phenomenon  which  can  only  be  de- 
scribed in  full  detail  by  the  time-dependent  scattering  theory.  This  implies  the  solution 
of  the  time-dependent  Schrodinger  equation  of  the  whole  system 


Hip  ( t ) = ih 


dip  ( t ) 

dt 


(2.17) 


when  the  Schroodinger  picture  is  adopted.  At  initial  time,  both  the  projectile  and  the 
target  are  described  by  the  wave  function  ip  (t  = — oo).  In  the  usual  case  of  a state - 
to-state  scattering  problem,  the  initial  wave  function  is  a stationary  eigenfunction  of 
the  whole  Hamiltonian  with  a trivial  time  dependence  through  a phase.  However,  it  is 
usual  in  some  type  of  simulations  not  to  employ  a single  state  but  some  linear  combina- 
tion of  state  forming  a wave  packet.  In  either  case,  at  a remote  final  time,  the  evolved 
wave  function  ip  (t  = oo)  can  be  decomposed  into  a linear  combination  of  the  station- 
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ary  eigenfunction  of  the  products.  Then,  the  coefficients  of  that  combination  ought  to 
have  a trivial  time-dependence  through  their  phase  again.  At  any  time  of  the  evolution 
-oo  <t<  oo,  all  the  properties  of  the  scattering  system,  including  final  time  properties 
such  as  cross  sections,  can  be  calculated  from  the  wave  function  ip  (t).  The  probability 
of  the  system  to  yield  by  reaction  some  type  of  products  at  a final  time  can  be  obtained 
by  projecting  the  final  wave  function  against  the  product  states.  As  shown,  the  scattering 
problem  also  demands  the  previous  solution  of  the  time-independent  eigenvalue  prob- 
lem of  both  reactants  and  products.  For  a time-independent  Hamiltionian,  the  totally 
formal  solution  of  this  problem  for  an  evolution  between  the  times  t\  < t2  is 


^(*2) 


u(h 

- ti)ip{ti) 

exp 

— ?:h  (t2  - ti) 

h 

^(<1) 


(2.18) 


where  the  central  unitary  operator  U (t2  — ti)  is  called  the  time-evolution  operator  or 
propagator. 

As  mentioned  the  most  detailed  description  of  a given  scattering  process  or  a chem- 
ical reaction  can  only  be  achieved  by  the  time-dependent  formulation  of  the  problem. 
However,  the  solution  of  a realistic  scattering  problem  is  far  more  difficult  to  find  in 
the  time-dependent  scheme  than  in  the  time-independent  one.  The  difficulties  posed 
by  the  time-dependent  methods  had  favored  the  adoption  of  the  time-independent  ap- 
proach until  the  1970s.  The  justification  of  how  essentially  time-dependent  phenomena 
such  as  scattering  processes  and  chemical  reactions  can  be  treated  time-independently 
is  rigorously  explained  in  the  more  advanced  literature  on  the  subject[3,  4],  Here,  it 
can  be  loosely  stated  that  the  eigenfunctions  ip  ( E ) of  the  time-independent  Schrodinger 
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equation  for  a scattering  problem 


H iP(E)  = Eip(E)  (2.19) 

contain  in  the  asymptotic  region  an  incoming  component  describing  the  reactants  and 
different  outgoing  components  describing  the  products.  Most  of  the  relevant  properties 
of  the  scattering  process  can  be  calculated  from  those  outgoing  components.  How- 
ever, all  the  intermediate  details  of  a scattering  process  or  a chemical  reaction  cannot  be 
obtained  in  this  way.  This  approach  is  an  extension  of  the  most  traditional  methods  to 
calculate  bound  states  in  isolated  molecules  to  the  case  of  several  colliding  molecules  in- 
volving some  unbound  translational  states.  This  connection  to  the  methods  in  electronic 
structure  theory  has  also  favored  the  bias  to  the  time-independent  schemes.  Because  of 
its  historical  precedence,  the  time-independent  methods  will  be  presented  first.  Further- 
more, the  time-independent  approach  directly  refers  to  the  stationary  description  of  re- 
actants and  products,  a problem  which  must  be  necessarily  addressed  before  attempting 
a time-dependent  solution.  It  is  also  important  to  remember  that  some  of  the  techniques 
further  developed  within  END  theory  have  their  origins  in  the  time-independent  theory. 

2.3.2  Fundamentals  of  the  Time-Independent  Theory 

For  the  simple  case  of  elastic  scattering  by  a localized  potential  V (r)  the  wave 
function  is  dependent  only  on  the  interparticle  distance  r.  After  Merzbacher[5]  and 
Sakurai[2]  motion  of  the  particle  is  described  by  the  Hamiltonian 

H = ^ + V = H0 + V 


(2.20) 
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where  H0  is  the  free  particle  Hamiltonian.  The  solution  to  the  time-independent  Schrodinger’s 
equation  is  then  given  by  the  Lippmann-Schwinger  equation[2] 

|^(±))  = \<t>)  + p 77 -|^(±)),  (2-21) 

£j  — 11  o — 26 


where  the  energy  term  in  the  denominator  has  been  made  complex  with  the  infinitesmal 
parameter  e in  order  to  ensure  continuity  and  the  ± indicates  outgoing  and  incoming 
solutions. 

Then  the  outgoing  solution  to  the  time-independent  Schrodinger  equation  has  the 
asymptotic  dependence 

</,(+)  (r)  r^°°  T/  + / (0)  (2.22) 

Here  '!'/  is  the  state  of  the  system  in  the  absence  of  the  scattering  potential  which  will  be 
determined  by  theoretical  or  experimental  considerations.  If  we  make  the  ansatz  that  the 
incident  wave  function  is  a plane  wave  which  represents  all  the  impact  parameters  with 
a given  momentum  kh  in  the  z direction  then  it  has  the  form  'I'/  = elkz . The  scattering 
amplitude  / (9)  governs  the  angular  distribution  of  the  scattering  which  is  outgoing 
due  to  the  factor  elkr . Other  forms  of  'I'/  may  be  applicable  in  certain  circumstances, 
for  example  if  the  dimension  of  the  incident  wave  function  were  small  compared  to  the 
effective  region  of  the  potential  requiring  a wave  packet.  Time-independent  wave  packet 
scattering  theory  is  fully  discussed  in  Goldberger[3]  and  Newton[6]. 

The  scattering  amplitude  is  related  to  the  DCS  through  its  definition  Eqn.  2.5  and 
the  evaluation  of  the  incident  and  scattered  flux  by 


dcTjj 

TlfT 


dH 


r2  ^ 


scattered 


l/(0)|dfi. 


(2.23) 


Jincident 


16 


Using  this  construction  the  problem  then  becomes  to  determine  / ( 9 ) for  a given  poten- 
tial function  V (r).  We  are  required  to  turn  our  attention  from  the  asymptotic  region 
to  the  range  where  the  potential  function  is  localized.  This  is  the  launching  point  for 
various  techniques  in  determining  the  form  of  scattering  amplitude. 

2.3.3  Scattering  Amplitude 

In  the  time  independent  formulation  given  by  Sakurai[2]  the  scattering  amplitude  is 

f(e)  = ~(2nf^(k'\VWM),  (2.24) 

— * 

where  k!  represents  the  propagation  vector  reaching  the  observation  point  0.  The  out- 
going wave  function  xp^  and  potential  V terms  in  Eqn.  2.24  are  the  fundamental 
unknowns  in  this  formulation.  We  will  leave  the  discussion  of  the  calculation  of  the 
potential  term  for  later.  The  determination  of  xp^  is  realized  through  several  different 
approximations  applicable  in  different  regimes.  Here  we  will  discuss  the  Bom  Approx- 
imation and  the  method  of  partial  waves  which  comprise  the  main  quantum  mechani- 
cal approximations,  then  we  will  address  the  semiclassical  JWKB  approximation,  and 
Eikonal  approximation. 

2.4  Semiclassical  Theory  of  Scattering 
2.4.1  Born  Approximation 

If  the  assumption  is  made  that  the  effect  of  the  scattering  potential  is  not  very  strong 
we  make  an  approximation  to  replace  xp^  with  the  free  particle  solution  (p.  In  this  way 
we  treat  the  effect  of  the  potential  V as  a perturbation  of  the  free  Hamiltonian.  In  the 
fashion  of  perturbation  expansion  we  can  calculate  the  fist-order  Bom  amplitude  which 
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is  denoted  by 


/(i)  (m  = -±-2™  [ ei*-k")  £'v  (f')  dV.  (2.25) 

J 47 r ft2  7 


We  can  generate  an  iterative  solutions  by  introducing  the  transition  operator  T such 
that 

V\rpM)  = T\4>)  (2.26) 

Using  the  Lippmann-Schwinger  equation  Eqn.  2.21,  we  obtain  an  expression  for  T 


T=V+V- 

E 


1 


H0  - it 


T 


(2.27) 


The  scattering  amplitude  then  can  be  written  as 

1 O m ->  -* 

f{9)  = -^-(2^w(k'\T\k).  (2.28) 


The  iterative  solution  for  T is  given  by 


T = V + V 


E — Ho  — it 


U-f  V 


E- Ho 


-V  E — Ho  — itV  + 


(2.29) 


it 


Correspondingly,  we  can  expand  / as  follows: 

OO 

/ (0)  = E /W  (0) , (2.30) 

71=1 

where  n indicates  the  number  of  times  U enters  the  expansion. 

2.4.2  Partial  Wave  Expansion 

In  this  method  we  assume  the  potential  V is  spherically  symmetric  which  allows 
us  to  expand  the  scattering  amplitude  in  terms  of  angular  momentum  eigenfunctions 
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and  the  phase.  As  a consequence  of  our  assumption  about  the  potential  the  transition 
operator  T (see  Eqn.  2.29)  commutes  with  L 2 and  L.  Therefore,  T is  a scalar  operator. 
Application  of  the  Wigner-Eckart  theorem  to  a scalar  operator  immediately  gives 

(E',  l m'\T\E,  l,  m)  = Tt  (E)  8w5mm. , (2.31) 


where  E is  the  initial  energy,  l and  m are  the  angular  momentum  eigenvalues.  Clearly 
T is  diagonal  in  both  l and  m.  Turning  our  attention  Eqn.  2.28  we  can  show 

1 9m  ->  _» 

/(«)  = 

= %'LY.  (*')  rr  (*)  ■ w) 

\K\  l m 

Since  we  have  assumed  the  incident  beam  is  moving  in  the  z direction  and  k is  outgoing 
wave  at  the  angle  9 we  see  that 


/ (0)  = £ (2*  + 1)  /i  ( E ) Pi  (cos  9)  , (2.33) 

i 


where/(  (E)  = -i tT[  ( E ) /\k\.  Considering  the  conservation  of  angular  momentum  and 
that  in  our  construction  the  flux  of  particles  is  conserved  it  is  possible  to  relate  the  / to 
the  phase  of  the  outgoing  wave 


elSi  sin  Si 


(2.34) 


Then  the  scattering  amplitude  is 


1 OO 

/W  = TnE(2,  + 1)eii'sin'5'flW 

\k\  i- o 


(2.35) 
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with  Si  real.  The  effect  of  this  expansion  is  to  change  the  problem  from  calculating  the 
scattering  amplitude  to  evaluation  of  the  effect  of  the  potential  on  the  phase  which  in 
some  cases  may  be  significantly  easier.  The  differential  cross  section  can  be  obtained 
just  as  in  Eqn.  2.23  and  the  total  cross  section  can  be  succinctly  given  by 

°tot  = J\f(e)\2<m 

Ajr  °° 

= 7T^(2/  + l)sin 2 Si.  (2.36) 

k2  i=o 

At  high  energies  we  may  make  an  additional  approximation  which  is  our  first  step 
into  the  realm  of  semiclassical  approximations.  That  additional  approximation  is  that  at 
high  energies  the  spacing  between  the  energy  eigenstates  becomes  quite  small.  In  this 

— a 

case  we  may  regard  the  angular  momentum  quantum  number  l = b\k\  as  continuous. 
We  can  then  take  the  sums  that  appear  in  Eqn.  2.36  and  Eqn.  2.35. 

2.4.3  JWKB 

Now  we  turn  our  attention  to  semiclassical  Jeffreys,  Wentzel,  Kramers  and  Bril- 
louin  (JWKB)  approximation.  The  wave  length  of  relative  motion  of  two  nuclei  in  a 
heavy-ion  reaction  is  small  and  semi-classical  methods  exploit  the  smallness  of  this 
quantity  to  obtain  simple  approximate  formulae  for  scattering  phase  shifts  and  partial 
wave  amplitudes.  Following  the  discussion  of  Brink[7],  the  JWKB  begins  by  recasting 
the  Schrodinger  equation  as 

+ k2  (x)  ■0  = 0,  (2.37) 

where  k 2 ( x ) = (2 m/h2^j  (E  — V (x)).  Hence,  we  can  interpret  k ( x ) as  the  local  wave 
number  and  p (x)  = hk  (x)  is  the  local  momentum  of  the  particle.  Substituting 


ip(x)  = elf^x)/h 


(2.38) 
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so  that  Eqn.  2.37  reduces  to 


V 


,2 


(2.39) 


The  JWKB  approximation  is  dependent  on  the  assumption  that  77  is  a slowly  varying 
function  of  x and  solving  equation  Eqn.  2.39  by  iteration.  The  second  iteration  gives 


2.4.4  Eikonal  Approximation 

The  Eikonal  Approximation^]  is  applied  to  cases  in  which  the  scattering  potential 
V varies  very  little  over  a distance  of  order  of  the  wavelength  A of  the  scatterer.  Note 
that  V itself  need  not  be  small  only  that  £>k;  which  is  in  contrast  to  the  applicability 
of  the  Bom  Approximation  (Pg.  16).  This  condition  allows  us  to  replace  the  exact  wave 
function  ipW  by  a semiclassical  wave  function 


where  S ( x ) can  be  interpreted  as  the  phase  of  the  wave  function.  Using  the  flux  we  can 
observe  that 


7)  (x)  ~ (x)  + -i  — 

h rC 


(2.40) 


and  therefore 


(2.41) 


^(+>  - eiS(*)/ft, 


(2.42) 


j ~ x/S 


(2.43) 
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allowing  us  to  interpret  \/S  as  the  “momentum”  of  the  wave  packet.  Using  Eqn.  2.42  in 
the  time  independent  Schrodinger  equation  leads  to  the  Hamilton-Jacobi  equation  for  S, 


(y£)2 

2 m 


+ V = E = 


h2k2 

2m 


(2.44) 


Here  again  we  have  transferred  the  difficult  problem  of  calculating  to  the  calcula- 
tion of  the  phase.  For  the  high  energy  regime  the  computation  of  S is  often  made  with 
the  additional  approximation  of  straight  line  trajectories  for  the  semiclassical  path  of  the 
scattered  wave  packet. 

2.5  Time-Dependent  Scattering  Theory 

While  the  time-independent  methods  have  dominated  the  development  of  scattering 
theory  for  most  of  the  history  of  quantum  mechanics  it  must  be  remembered  that  scat- 
tering is  a time  dependent  process.  Time-dependent  theories  were  not  studied  in  earnest 
until  the  advent  of  modem  computers  in  the  1970s  . One  of  the  main  reasons  being  that 
the  time-dependent  equations  are  posed  far  more  demanding  then  the  time-independent 
ones.  This  young  branch  of  scattering  theory  is  not  so  developed  and  established  as 
its  counterpart  making  a concise  summary  more  difficult.  Several  reviews  have  been 
compiled  among  them  are  Deumens  et  al . [8]  and  Kroger[9]. 

The  time  dependent  methods  can  be  roughly  organized  into  two  camps.  The  first 
being  the  case  where  the  time  dependence  is  reserved  only  to  the  nuclear  degrees  of 
freedom.  In  the  second  one,  the  time  dependence  is  given  to  all  the  degrees  of  free- 
dom, nuclear  and  electronic.  The  first  approach  is  time-dependent  nuclear  dynamics 
on  one  or  more  predetermined  potential  energy  surfaces  (PES).  In  calculating  the  PES 
the  electronic  degrees  of  freedom  have  been  included  before  the  dynamical  equations 
are  solved  and  propagated.  In  the  second,  the  time  evolution  of  electrons  and  nuclei  is 
performed  simultaneously.  However,  this  calculation  is  very  intensive  if  PES  are  not 
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used  so  practical  methods  have  so  far  been  mostly  restricted  to  classical,  semiclassical, 
or  quasiclassical  descriptions  of  the  nuclear  degrees  of  freedom. 

2.5.1  Wave-Packet  Propagation:  Exact  Methods 

Methods  for  which  the  numerical  integration  of  the  nuclear  time-dependent  Schrodinger 
equation  is  performed  directly  are  called  exact[10].  Although  these  methods  contain  ap- 
proximations for  the  structure  of  the  Hamiltonian  and  other  numerical  approximations 
they  can  be  calculated  to  any  desired  level  of  accuracy.  The  scattering  is  simulated  and 
the  use  of  a wave  packet  either  represented  in  terms  of  discretized  grid  of  points  or 
expanded  in  a convenient  basis  set. 

2.5.2  Approximate  time-dependent  methods 

Time-dependent  self  consistent  field.  For  these  methods  the  exact  time  evolution 
is  replaced  with  numerically  easier  approximations  and  models.  One  such  method  is 
the  Time-Dependent  Self  Consistent  Field  (TDSCF)  method  introduced  by  Bisseling  et 
al.[l  1],  The  nuclei  are  represented  by  Hartree  products  of  wave  packets  on  a grid.  Each 
nucleus  is  then  moving  in  the  average  field  of  the  others.  Further  developments  of  this 
method  add  several  exit  channels  and  include  some  correlation  effects  for  the  nuclear 
dynamics  by  adding  a multi-reference  scheme  of  Hartree  products.  These  methods  can 
be  used  for  scattering,  photodissociation,  atomic  transfer  and  vibrational  dissociation. 

Trajectory  Surface  Hopping.  The  trajectory  surface  hopping  method  (TSH)  uses 
a probabilistic  description [12,  13]  to  make  the  solution  of  the  time-dependent  equations 
easier.  Instead  of  solving  the  coupled  set  of  equations  at  all  nuclear  configurations,  a 
hopping  probability  for  the  electronic  configuration  is  computed  at  each  instant.  Thus 
the  equation  for  the  classical  nuclei  depends  only  on  the  hopping  probability  between 
the  states  involved,  usually  just  two.  Direct  calculation  of  the  dynamics  of  the  electronic 
states  and  phase  is  ignored  in  the  TSH  method[14]. 
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Car-Parrinello.  A very  popular  method  for  dynamics  is  the  the  Car-Parrinello  (CP) 
method[15].  The  original  method  was  an  algorithm  to  solve  minimization  problems 
which  lead  the  the  solution  of  eigenvalue  problems.  While  this  method  could  be  used 
to  solve  any  number  of  eigenvalue  equations  it  has  primarily  been  applied  using  the 
equations  derived  from  Density  Functional  Theory  (DFT) 

Time-dependent  Hartree-Fock.  This  large  class  of  methods  is  based  on  the  Hartree- 
Fock  equation  for  the  dynamics[16]  and  utilizes  semiclassical  or  classical  descriptions 
for  the  atomic  nuclei.  These  methods  consider  an  explicit  dynamical  description  of  the 
electronic  state.  Sometimes  the  full  ab  initio  Flamiltonian  is  considered,  in  other  cases 
model  Hamiltonians  are  set  up  to  drive  the  dynamics.  The  coupling  between  the  elec- 
trons and  nuclei  in  these  models  is  through  the  average  potential  energy  surface.  The 
nuclei  and  the  electrons  affect  each  other  only  through  their  instantaneous  positions  in 
the  Fock  operator.  The  result  of  this  that  high  energy,  non-chemical  scattering  is  not 
properly  handled.  The  ad  hoc  solution  for  this  problem  is  the  addition  of  electron  trans- 
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CHAPTER  3 

END  THEORY  OF  TIME  DEPENDENT  MOLECULAR  DYNAMICS 


Electron  Nuclear  Dynamics  (END)[17]  is  becoming  a complete  general  theory  of  molec- 
ular processes.  It  is  one  theory  of  a few  among  the  dynamics  community  in  that  its 
foundations  are  based  in  the  Time  Dependent  Variational  Principle  (TDVP).  This  novel 
approach  generates  a theoretical  framework  to  calculate  the  full  dynamics  of  molecu- 
lar collisions  and  dynamics  without  the  use  of  the  Bom-Oppenheimer  approximation 
thus  allowing  the  study  of  nonadiabatic  interactions.  The  richness  of  the  theory  can  be 
implemented  in  a hierarchy  of  complexity  ranging  from  the  semiclassical  to  the  fully 
quantum.  At  all  levels  of  the  realization  of  this  theory  there  is  a fertile  ground  of  inter- 
esting physical  and  chemical  problems  to  investigate. 

Currently  the  END  theory  is  implemented  in  the  ENDyne[18]  code.  The  level  of 
implementation  uses  classical  nuclei  in  the  narrow  wave  packet  limit  and  a single  com- 
plex determinant  to  describe  the  electrons.  Even  at  this  basic  level  of  approximation 
several  applications  have  shown  that  END  can  capture  much  of  the  physics  of  nonadi- 
abatic molecular  processes[8,  19,  20,  21,  22,  23,  24,  25].  It  has  also  become  clear  that 
there  are  limitations  to  this  level  of  the  theory.  Further  developments  in  implementation 
are  in  progress. 

3.1  The  END  Theory 

The  starting  point  for  the  END  theory  is  the  TDV[26]  which  is  the  quantum  me- 
chanical analogue  of  the  classical  principle  of  least  action  sometimes  called  “Hamil- 
ton’s Principle”.  The  application  of  the  principle  of  least  action  to  a physical  problem 
often  leads  to  more  computationally  tractable  equations  of  motion  than  for  example  the 
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Schrodinger’s  equation.  Here  we  will  introduce  the  quantum  mechanical  action  A and 
demonstrate  its  connection  to  the  Schrodingers  equation  by  introducing  a general  set 
of  complex  variational  parameters.  We  will  then  derive  the  equations  of  motion  using 
the  overlap  matrix  and  the  energy.  Finally  we  will  discuss  the  implementations  of  END 
specifically  the  implementation  featured  in  the  ENDyne  code. 

The  quantum  mechanical  action  is  given  by 

A = [t2  Ldt,  (3.1) 

Jti 


with  the  quantum  mechanical  Lagrangian  defined  as  ( h = 1) 


(3.2) 


Here  H is  the  quantum  mechanical  Hamiltonian  of  the  whole  system.  Note  the  “left 
acting  derivative”  used  in  the  Lagrangian.  It  is  defined  for  a differentiable  function  / 
as 


fl  = If. 

J dt  dV 


(3.3) 


This  definition  is  necessary  if  we  want  to  promote  Eqn.  3.2  to  an  operator  in  Hilbert 
space,  L — » 0,  and  ensure  it  is  Hermitian.  We  make  the  ansatz  that  4/,  the  wave 
function  for  the  whole  system,  depends  on  a set  of  complex  variational  parameters 
£ = {Ci;  (2,  (3,  ■ ■ ■ 1 Cm}  which  depend  formally  on  time,  Q = £(t).  We  introduce 
the  column  vector  |£)  as  the  representation  of  T = 'F(C)  = |C)  in  terms  of  the  complex 
parameters  {£i}-  Here  we  observe  that  these  complex  parameters  along  with  their  time 
derivatives  serve  as  generalized  positions  and  velocities  for  our  quantum  mechanical 
Lagrangian.  Thus 


£(*,**)  = L(C,C*,C,a, 


(3.4) 
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where  £ = ^ are  the  generalized  velocities.  The  time  evolution  of  the  system  is  deter- 
mined using  the  TDVP  which  requires  that  the  action  be  an  extremum  with  respect  to 
its  parameters  so 


SA  = 5 


rh 

/ Ldt  = 0, 

J 


(3.5) 


subject  to  the  boundary  conditions 


8\V)  = <5(\&|  = 0 


(3.6) 


at  t = ti  and  t2[21].  Substituting  our  expressions  and  employing  the  above  boundary 
conditions  and  integration  by  parts  gives  us 


6A 


5 [hdt  (c|0|C)(ClC)  = o 

Jti 


mm  (cieio 

(CIO  (CIO2 


(3.7) 


Since  'P,  and  equivalently  (,  can  be  varied  throughout  the  full  Hilbert  space  the  integrand 
of  Eqn.  3.7  must  satisfy 

[4  - = ®o  <3-8> 

If  we  explicitly  consider  the  phase,  7,  of  the  wave  function  by  writing  |0  — > e*7|0  we 
can  derive  the  following  expression  for  the  right  hand  side  of  Eqn.  3.8 


(C|e-J7©ei7|0  = 0, 


(C|e  7)  - i(-i i))en\()  + (CI@|C)  = 0, 


(CIO7  + (Cieio  = 0, 


(3.9) 
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Substituting  these  results  in  Eqn.  3.8  gives 


mh  -h]  |c)  = 7lC>. 
Ht-H]  10  - OIO  = o 


=*  [iftg-ff]  ei7|C)  = 0. 


(3.10) 


This  result  is  the  Schrodinger  equation.  The  preceding  derivation  provides  the  justi- 
fication for  the  form  of  the  action  presented  in  Eqn.  3.4.  If  we  choose  to  specify  a 
constrained  form  for  the  wave  functions  'k  then  we  would  restrict  ourselves  to  particular 
regions  of  the  Hilbert  space  resulting  in  a dynamical  evolution  that  would  be  an  approx- 
imate solution  to  the  Schodinger  equation.  This  idea  forms  the  kernel  from  which  all 
realizations  of  the  END  theory  are  derived.  In  choosing  the  set  of  parameters  ( and  the 
form  of  the  Hamiltionian,  H,  we  are  able  to  derive  equations  of  motion  that  are  relevant 
to  many  different  molecular  processes.  The  primary  advantage  of  this  theory  is  in  the 
fact  that  we  have  not  constrained  the  dynamics  of  the  system  except  in  our  choice  of 
parameterization.  Hence  we  are  not  limited  to  investigating  only  adiabatic  processes  as 
molecular  studies  using  PES  often  are. 

Even  before  making  any  assumptions  about  the  form  of  we  can  derive  a set  of 
general  equations  of  motion.  Let  us  define  the  S = S( (,(*)  = (CIO  and  the  energy 
E = E{ C,  C*)  = (C|H|C)/(C|C)- Then  Eqn.  3.7  can  be  expressed  in  ( representation  as 


(3.11) 
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Since  S(p  and  6Q  are  independent  variations  then  two  sums  in  the  integrand  must  be 
separately  zero  resulting  in 


P 


dCa 


^ = it' 


(3.12) 


where  we  define  the  complex  Hermitian  matrix  C = {CQ/g}  = d2  In  S/d(*d(p.  The 
general  END  dynamical  equations  are  given  by  Eqn.  3.12  and  they  may  be  recast  in 
matrix  form  as 

( n n \ ^ ^ ^ M 


V 


c 0 
0 -c 

We  define  a generalized  Poisson  bracket  by 


dE 

V 9Ca  ) 


(3.13) 


UiS}  = -*£ 

a,P 


d 'f_c _i  0g  _ dg  x df 

d(a  a/3dQ  d(a  a^dQ 


(3.14) 


with  / (£,  £*)  and  g ((,  (*)  being  differentiable  functions.  It  follows  that  ( = {C,  E}  and 
(*  — {£*,  E}.  So  we  can  state  that  the  generalized  Poisson  brackets  show  that  the  time 
evolution  of  the  wave  function  parameters  are  governed  by  Hamilton-like  equations. 
Moreover,  the  additional  relations 


{C>,  Oj}  = {C  Q}=  0 and  {C„,  q}  = (C  0}  = -iC~l  (3.15) 

show  that  £*  and  ( behave  as  “classical”  coordinates.  Observe  that  if  the  matrix  C were 
the  unit  matrix  then  Eqn.  3.15  would  reduce  to  the  commutation  rules  for  a flat  phase 
space.  Our  generalization  has  lead  to  a curved  phase  space. 


29 


3.2  The  END  Wave  function 

At  the  time  of  this  writing  only  the  first  nontrivial  level  of  END  theory  has  been  im- 
plemented in  the  code  ENDyne[18],  In  this  model  the  unnormalized  END  wave  function 
in  the  lab  frame  is  written  as 


^ END  (X,  X,  t)  Fnud  (X;  R(t),  P(tj)  fei  (x;  z (t),  R(t))  x 


exp 


hl total  (t) 


(3.16) 


where  X and  x are  the  nuclear  and  electronic  coordinates,  R(t),  P(t)  and  z (t)  are  the 
variational  parameters,  and  7 total  (t)  is  the  total  phase.  The  time-dependence  of  the  wave 
function  comes  only  through  the  parameters  and  the  phase.  We  describe  this  wave  func- 
tion as  the  quasiclassical  single  determinental  END  wave  function  (QCSD  END)  and 
it  has  three  important  parts:  the  nuclear  wave  function  Fnuci  (X;  R(t),  P(£)),  the  elec- 
tronic wave  function  fet  (x;  z (t),  R(t)^,  and  the  phase. 

The  nuclear  part  consists  of  a simple  product  of  unnormalized,  so-called  “frozen” 
Gaussian  wave  packets  $ (Xk ; Rk  (t) , Pk  (£)) 

N Till  cl 

Fnucl  (X;  R(t),  P(t))  = A ^ (X,;  Rk  (t) , Pk  (t))  , (3.17) 

k= 1 


where 


(3.18) 


The  variational  parameters  R(t),  P(t),  and  the  constant  ak  are  assumed  to  be  real.  The 
construction  of  these  wave  functions  establishes  the  connection  of  these  parameters  with 
average  position  and  momenta  of  the  wave  packets.  The  constant  ak  can  be  related  to 
standard  deviation  or  width  of  the  Gaussian  functions. 
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The  physical  meaning  of  the  previous  parameters  and  constants  imply  that  each  nu- 
clear wave  packet  can  be  seen  as  parametric  plane  wave  function 


Furthermore,  since  the  constants  ak  are  time-independent  and  the  evolution  is  through 
the  TDVP,  each  nuclear  Gaussian  will  remain  a Gaussian  during  the  evolution  without 
changing  shape.  The  nuclear  wave  packets  only  change  in  their  positions  and  momenta. 
Contrasting  this  method  with  unconstrained  propagation  of  nuclear  wavepackets,  we 
observe  the  QCSD  END  nuclear  wave  functions  is  not  allowed  to  spread,  change  shape 
or  split  into  seperate  packets.  The  frozen  wavepackets  are  limited  in  this  regard  but  have 
great  computational  advantage  over  the  general  method.  The  use  of  frozen  Gaussian 
wavepackets  has  been  widespread[3]  and  is  considered  quite  valid  if  one  chooses  to 
ignore  tunneling.  In  most  molecular  experiments  the  scattering  process  happens  is  such 
that  the  width  of  the  nuclear  wavefunction  can  be  considered  negligable.  Additionally, 
this  level  of  wave  function  for  the  nuclei  allows  a description  of  the  internal  motion  of 
the  molecular  components  of  the  scattering  process. 

The  electronic  function  fei  fx;  z (£),  R(t)^j  is  a single-determinant  wavefunction  through- 
out the  time  evolution.  This  type  of  function  has  been  long  established  in  the  TDHF 
theory  of  nuclear  many-body  theory  [28,  29,  30].  Electron  Nuclear  Dynamics  theory  has 
an  attractive  feature  in  its  ability  to  treat  both  the  nuclear  and  electronic  degrees  of  free- 
dom simultaneously.  Additionally,  this  single  determinant  is  related  to  either  restricted 
or  unrestricted  Hartree-Fock  description  of  the  reagents  and  products  ground  states.  The 


exp  j^Pk  ( t ) • Xfc  - Rk  (t) 


(3.19) 


centered  around  the  position  R(t)  by  a simple  Gaussian  factor 


(3.20) 
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difficult  problem  of  constructing  a wavefunction  which  remains  as  a single  determinant 
while  being  evolved  was  solved  by  Thouless[29,  30],  The  electronic  QCSD  END  wave 
function  is  presented  following  the  original  derivation  by  Thouless. 

It  is  meaningful  to  divide  the  set  of  K spin  orbitals  into  N occupied  and  K — N 
unoccupied  spin  orbitals.  The  linear  array  q refers  to  the  set  of  all  all  spin  orbitals.  Then 
q*  and  q°  refer  to  the  occupied  and  unoccupied  orbitals,  respectivly.  The  atomic  spin 
orbitals  may  also  be  partitioned  into  two  sets.  We  write  them  as 

<f>  = W <t>°)  ■ (3.21) 


Similarly,  for  the  matrices  in  the  equations  four  sub-blocks  are  identified:  the  occupied 
N x N and  unoccupied  (. K — N)  x (K  — N)  diagonal  blocks;  and  the  upper  and  lower 
off  diagonal  blocks.  The  transformation  to  the  molecular  basis  is  then 


/ 

(r  r)  = (0*  <n 

v 


w* 

kEv 


W>  ^ 

w° 


(3.22) 


The  right-angle  > denotes  the  upper  off-diagonal  block,  and  the  down-angle  V denotes 
the  lower  off-diagonal  block  which  has  a vertical  rectangular  shape.  The  purposes  of 
this  notation  is  to  avoid  many  messy  index  manipulations.  As  a mnemonic  we  reserve 
p,  q,  r for  indices  running  over  the  particle  or  unoccuupied  range  N + 1,  . . . , K\  and 
h,  g,  f are  for  the  hole  or  occupied  range  1,  . . . , N.  Both  the  atomic  <j>  and  molecular 
0 depend  parametrically  on  the  Gaussian  parameter  R ( t ). 
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Introducing  a general  unitary  transformation  U of  the  orthonormal  molecular  spin 
orbital  basis  one  can  write 


(c#t  cot)  = (Vf  hot) 


1 U'  U>  ' 


[7V  U° 


(3.23) 


for  the  basis  field  operators,  and  obtain  for  a determinental  state[31] 

|4'>  = |c-t)  =n£ic;'|»ac)  =anfl+  E E K'u*u^'b'i! ) ft  V I «“> 

h=l  \ p=N+ 1 k= 1 / /=1 

where  the  invariance  of  a determinantal  wave  function  under  a linear  transformation 
of  its  occupied  spin  orbitals  up  to  a constant  a has  been  used.  Defining  the  Thouless 
parameters  zpq,  EjtLi  = zph,  and  the  reference  state 


N 


= 


6#t\  = n^f  \va°) 


(3.24) 


i=i 


one  can  write  the  unnormalized  Thouless  determinant 


Z,  R)  = fei  (x;  z (t),  fl(t))  as 


Z ,R)=  nLi  (l  + Ef=jv+i  tfzphbh)  I ^0  ) 

= exp  (E/ll  Ep=n+ 1 Zphbpbl)  | 0 ) . 


(3.25) 


Here  the  nilpotence  of  the  operators  b^ 6*  has  been  used.  The  parameters  zpq  are  com- 
plex numbers.  The  determinantal  wavefunction  of  this  state  is  det  Xh  (fn)  in  terms  of 
the  occupied  dynamical  orbitals 

K 

Xh='lph+  Pzph ■ 

p=./V-fl 


(3.26) 


33 


Z,  R)  = fei  (x;  z (t),  R(tj)  is  the 


These  orbitals  are  not  orthogonal.  The  function 
single-determinant  wave  function  for  the  QCSC  END  theory.  If  the  arbitrary  reference 
state  | \&o ) is  selected  as  the  HF  ground  state  of  the  system  then  the  Thouless  determinant 
contains  this  state  plus  all  its  excitations.  Finally,  the  QCSD  END  total  phase  7 totai  (t) 
is  the  quantum  action  corresponding  to  the  QCSD  END  wave  function  'end  (X,  x,  t). 


CHAPTER  4 

OVERVIEW  OF  ROVIBRATIONAL  REACTION  DYNAMICS 


In  this  chapter  we  provide  a brief  overview  of  the  classical  and  quantum  mechanical 
rovibrational  dyanamics.  First,  we  review  classical  vibrations.  Second,  we  discuss 
the  quantum  mechanical  description  of  the  harmonic  oscillator.  The  classical  theory 
of  rotations  is  examined.  Finally,  the  quantum  mechanical  description  of  rotations  is 
presented. 


A system  is  said  to  be  in  equilibrium  when  the  generalized  forces  acting  on  the 
system  vanish 


The  potential  energy  therefore  has  an  extremum  at  the  equilibrium  configuration  of  the 
system,  §01,  q02,  . . . , qon • If  the  initial  configuration  is  at  the  equilibrium  position  with 
zero  initial  velocities  q^,  then  the  system  will  remain  in  equilibrium,  i.e.  like  a pendulum 
at  rest.  In  the  case  of  a stable  equilibrium,  one  for  which  small  displacements  result  in 
small  bounded  motion  about  the  resting  position,  all  functions  can  be  expanded  in  terms 
of  a Taylor  series  and  we  retain  only  the  lowest  order  terms.  We  naturally  cast  our 
system  in  terms  of  deviations  r]l  from  equilibrium: 


4.1  Classical  Theory  of  Vibrations 


(4.1) 


q%  qoi 


(4.2) 
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Using  r?i  as  the  new  generalized  coordinates  of  the  motion  we  expand,  the  potential 
energy  about  q0l  we  obtain 


V (qu  q2,  . . . , qn)  = V {q0 1,  q0 2,  ■ ■ • , qon)  + 


d2V  \ 
dqidqj  J 


ViVj  + ■ ■ ■ 
o 

(4.3) 


Where  we  have  imposed  the  Einstein  summation  convention.  As  we  see  from  Eqn.  4.1 
the  linear  term  vanishes  and  we  shift  the  zero  of  potential  energy  to  eliminate  the  first 
term  in  the  series.  We  are  then  left  with  the  quadratic  terms  as  the  first  approximation  to 


/dU\  1 l d2V 
\dqldqj 


ViVj  = ^VrjViVj, 
0 z 


(4.4) 


where  the  second  derivatives  of  V have  been  designated  by  the  constants  V^. 

A series  expansion  can  also  be  obtained  for  the  kinetic  energy.  The  kinetic  energy 
is  a homogeneous  quadratic  function  of  the  velocities  when  the  generalized  coordinates 
do  not  explicitly  depend  on  time: 


-f  2 ^ijQiQj  2 j Vi  Vj  ■ 


(4.5) 


The  coefficients  mtj  are  functions  of  the  coordinates  qk,  but  they  may  be  expanded  in  a 
Taylor  series  about  the  equilibrium  configuration  just  as  we  did  for  V . Keeping  only  the 
quadratic  terms,  we  can  therefore  write  the  kinetic  energy  as 


T = 2TbW 


(4.6) 
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We  should  observe  that  both  Vtj  and  Ttj  are  symmetric  since  individual  terms  are  unaf- 
fected by  interchange  of  indices.  The  Lagrangian  is  given  now  given  by 


Using  the  Vi  as  the  generalized  coordinates,  the  Lagrangian  leads  to  the  following  n 
equations  of  the  motion: 


where  explicit  use  has  been  made  of  the  symmetry  property  of  the  VtJ  and  coeffi- 
cients. We  then  must  solve  this  set  of  simultaneous  differential  equations  involving  all 
of  the  coordinates  Vi  to  obtain  the  motion  near  the  equilibrium. 

The  equations  of  motion  are  linear  differential  equations  with  constant  coefficients 
leading  us  to  make  the  ansatz  that  the  solution  has  the  form 


Here  Ca*  gives  the  complex  amplitude  of  the  oscillation  for  each  coordinate  Vi,  the 
factor  C being  introduced  for  convenience  as  a scale  factor.  Clearly,  it  is  the  real  part  of 
Eqn.  4.9  that  corresponds  to  the  real  motion.  Substituting  our  ansatz  into  the  equations 
of  motion  leads  to 


(4.7) 


TijVj  + VijVj  — 0) 


(4.8) 


Vi  = Ca,e  lbJt 


(4.9) 


(4.10) 
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These  equations  constitute  n linear  homogeneous  equations  for  the  a*’s  and  conse- 
quently can  have  a solution  only  if  the  determinant  of  the  coefficients  vanishes: 


Vn-u2Tn  Vu  - cu2Tu 

Vn  ~ U2T2i  V12  — U2T22 

Vn-u2T2l 


(4.11) 


We  have  in  effect  an  algebraic  equation  of  the  nth  degree  for  u2,  and  the  roots  of  the 
determinant  provide  the  frequencies  for  which  Eqn.  4.9  represents  a correct  solution 
to  the  equations  of  the  motion.  Each  of  these  values  of  cu2  the  equations  maybe  solved 
for  the  amplitudes  of  a*,  or  more  precisely,  for  n — 1 of  the  amplitudes  in  terms  of  the 
remaining  a;. 

The  general  solution  of  the  equations  of  motion  may  now  be  written  as  a summation 
over  an  index  k: 

Vl  = Ckalke~'“kt,  (4.12) 

there  being  a complex  scale  factor  Ck  for  each  resonant  frequency.  We  should  note 
that  we  are  only  interested  in  the  real  part  of  our  equations  of  motion,  consequently 
we  ignore  the  negative  values  of  u>.  It  is  possible  to  transform  the  t?,  to  a new  sort  of 
generalized  coordinates  that  are  all  simple  periodic  functions  of  time-a  set  of  variables 
known  as  the  normal  coordinates. 

We  define  a new  set  of  coordinates  Q related  to  the  original  coordinates  77*  by  the 
equations 

Tji  (1/j  , (4.13) 

or  in  linear  algebra  terms 


V = M- 


(4.14) 
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Importantly,  A diagonalizes  V by  a congruence  transformation  and  the  potential  there- 
fore reduces  simply  to 


(4.15) 


Similarly,  the  kinetic  energy  transforms  to 


(4.16) 


Returning  to  the  Lagrangian,  we  know  have 

L = l (CiCi  - W,2C,2) 


(4.17) 


so  that  the  Lagrange  equations  for  Q are 


Ci  + ^iCi  — 0. 


(4.18) 


The  immediate  solutions  are 

C i = Cie~iUit.  (4.19) 

Each  of  the  new  coordinates  is  thus  a simple  periodic  function  involving  only  one  of  the 
resonant  frequencies.  So  it  is  therefore  customary  to  call  the  C’s  the  normal  coordinates 
of  the  system. 

Individually  each  normal  coordinate  corresponds  to  a vibration  of  the  system  with 
only  one  frequency,  and  these  component  oscillations  are  called  the  normal  modes  of 
vibration.  All  of  the  particles  in  each  mode  vibrate  with  the  same  frequency  and  with 
the  same  phase;  the  relative  amplitudes  being  determined  by  the  matrix  elements  aik. 
The  complete  motion  is  then  built  up  out  the  the  sum  of  the  normal  modes  weighted 
with  appropriate  amplitude  and  phase  factors  contained  in  the  Ck  s. 
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Harmonics  of  the  fundamental  frequencies  are  absent  in  the  complete  motion  essen- 
tially because  of  the  stipulation  that  the  amplitude  of  oscillation  be  small.  We  are  then 
allowed  to  represent  the  potential  as  a quadratic  form,  which  is  characteristic  of  simple 
harmonic  motion.  The  normal  coordinate  transformation  emphasizes  this  point,  for  the 
Lagrangian  in  the  normal  coordinates  is  seen  to  be  the  sum  of  Lagrangians  for  harmonic 
oscillators  of  frequencies  ujk.  We  can  thus  consider  the  complete  motion  for  small  os- 
cillations as  being  obtained  by  exciting  the  various  harmonic  oscillators  with  different 
intensities  and  phases. 

4.2  Quantum  Mechanical  Theory  of  Vibrations 

The  quantum  mechanical  harmonic  oscillator  is  a superb  pedagogical  tool  as  it  is  a 
system  that  can  be  exactly  solved  in  classical  and  quantum  theory  but  it  is  also  a system 
of  enormous  physical  relevance.  Any  system  displaced  by  small  amounts  near  a config- 
uration of  stable  equilibrium  may  be  described  either  by  an  oscillator  or  by  a collection 
of  decoupled  harmonic  oscillators.  The  dynamics  of  a collection  of  noninteracting  os- 
cillators is  no  more  complicated  than  that  of  a single  oscillator  aside  from  the  ./V-fold 
increase  in  degrees  of  freedom.  While  addressing  the  solution  of  the  oscillator  we  are 
actually  confronting  the  general  problem  of  small  oscillations  near  equilibrium  of  an 
arbitrary  system. 

4.2.1  The  Quantum  Mechanical  Harmonic  Oscillator 

The  canonical  example  of  a single  harmonic  oscillator  is  a mass  m coupled  to  a 
spring  of  force  constant  k.  Small  deviations  in  the  position  x will  exert  the  force  given 
by  Hooke’s  law,  F = ~kx,  and  produce  a potential  V = \kx2.  The  Hamiltionian  for 
this  system  is 

H = T + Vr  = — + \muj2x2 
2m  2 


(4.20) 
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where  oj  = (k/m)  2 is  the  classical  frequency  of  oscillation.  Any  Hamiltonian  that  is 


quadratic  in  the  coordinates  and  momenta  will  be  called  the  harmonic  oscillator  Hamil- 
tonian. The  mass-spring  system  is  just  one  among  the  following  family  of  systems  de- 
scribed by  the  oscillator  Hamiltonian.  A particle  moving  in  a potential  V (x)  is  placed 
at  one  of  its  minima  x0,  it  will  remain  there  in  a state  of  stable,  static  equilibrium.  Con- 
sider the  dynamics  of  this  particle  as  it  fluctuates  by  small  amounts  near  x = x0.  The 
potential  may  be  expanded  as  a Taylor  series  as  in  the  classical  model.  For  small  oscil- 
lations, we  may  again  neglect  all  but  the  leading  term  and  arrive  at  the  potential  Eqn. 
4.20.  Therefore,  we  identify  d2V/dx2  with  k = mo;2.  A system  of  harmonic  oscilla- 
tors can  be  formed  by  simple  addition  of  the  individual  Hamiltonians  and  decoupling 
them  through  a coordinate  transformation.  In  general,  a system  of  N coupled  harmonic 
oscillators  can  be  decoupled  into  N separate  harmonic  oscillators. 

4.2.2  Quantization  of  the  Harmonic  Oscillator 

We  focus  our  attention  on  the  time-independent  Schrodinger  equation  eigenvalue 
problem: 


A clever  method  due  to  Dirac  allows  us  to  work  in  the  energy  basis  without  having  to 
know  ahead  of  time  the  operators  X and  P.  From  the  basis  independent  commutation 
relation 


(4.21) 


[x,P  ] = ih, 


(4.22) 


we  are  motivated  to  introduce  the  operator 


(4.23) 
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and  its  adjoint 


+ / muj\  2 . / 1 \ 2 

a = \ 2h)  x~l  Gwi)  p 


They  satisfy  the  commutation  relation 


(4.24) 


a,  af  = 1, 


(4.25) 


and  is  related  to  the  harmonic  oscillator  Hamiltionian  so  that 


H = f a^a  + - ) hu>. 


(4.26) 


For  simplification  let  us  define  H = H /hu  and  e = E/hu.  Two  important  relations  can 
be  immediately  obtained: 


a,  H 
af , H 


(4.27) 


The  utility  of  a and  at  is  that  given  a particular  eigenstate  they  can  be  used  to  generate 
others.  For  example 


Ha  | e ) 


(aH 

(aH 

(c- 


a,H])|e) 
a)  |e) 


l)a|e) 


(4.28) 


We  can  conclude  that  a | e ) is  an  eigenstate  of  H with  eigenvalue  e — 1 thus 


a | e } = Ct  \ e — 1 ) 


(4.29) 
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where  Ce  is  a constant  and  | e ) and  | e — 1 ) are  normalized  eigenkets. 
Similarly  we  observe  that 


Haf  | e ) 


(a*H  + a^  | e) 


(e  + 1)  af  | e) 


(4.30) 


so  that 

af|e)  = Ce+1|c  + l).  (4.31) 

Clearly  this  is  why  one  refers  to  a and  a^  as  creation  and  annihilation  operators.  We  are 
led  to  conclude  that  if  e is  an  eigenvalue  of  H,  so  are  e±l,  e±2,  ..,  e±oo.  However,  we 
must  remember  that  the  quadratic  terms  in  H make  it  a positive  operator  with  positive 
eigenvalues.  Consequently,  there  must  exist  a state  | e0  ) that  cannot  be  lowered  further: 

a | e0  ) = 0 (4.32) 

and 

H|e0)  = ^.  (4.33) 

It  is  convenient  to  define 

en  = (n  + |)  , n = 0,  1,  2,  . . . 

En  — (n  + n = 0,  1,  2,  . . . (4.34) 

Since  there  are  no  degeneracies  in  one  dimension  these  constitute  all  of  the  eigenstates 


of  H. 
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With  some  manipulation  it  can  be  shown  that  the  coefficients  are  given  by 


(4.35) 


where  </>  is  an  arbitrary  phase.  The  energy  eigenstates  can  be  succinctly  written  in  terms 
of  the  creation  operator 


Returning  to  the  position  basis  with  X — > x and  P — > d/dx  we  can  observe  that  the 
inversion  of  Eqn.  4.23  and  Eqn.  4.24  leads  to 


4.3.1  Rotations 

When  we  discuss  rotations  we  do  so  in  the  context  of  a rigid  body.  A rigid  body  is 
a system  of  mass  points  subject  to  the  constraints  that  the  distances  between  all  pairs 
of  points  remain  constant  throughout  the  motion.  This  is  something  of  an  idealization, 
especially  in  the  case  of  molecular  motion,  but  it  allows  us  to  discuss  the  important 
aspects  of  rotational  kinematics  and  dynamics.  In  the  case  of  a rigid  body  with  N 
particles  it  can  at  most  have  3 N degrees  of  freedom.  Although,  as  we  see  by  definition 
there  exists  a set  of  constraints.  These  constraints  serve  to  reduce  the  number  of  degrees 
of  freedom  greatly.  We  only  need  to  establish  the  position  of  just  three  non-collinear 


(4.36) 


n ) = lf>n  (z) 


(4.37) 


where  Hn  ( y ) are  the  Hermite  polynomials. 


4.3  Classical  Theory  of  Rotations 
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particles  of  the  body  to  define  its  location.  This  gives  us  nine  degrees  of  freedom, 
except  we  are  not  free  to  alter  the  distances  between  the  particles  which  reduces  the 
total  number  to  just  6 degrees  of  freedom.  Three  of  which  describe  the  location  of  the 
rigid  body  and  the  other  three  describe  the  rotations. 

Rotations  of  rigid  bodies  are  thought  of  as  orthogonal  transformations  where  the 
coordinates  of  the  position  of  the  particles  in  the  rigid  body  are  transformed  into  another 
set  of  coordinates.  The  transformation  can  be  written  as: 

3 

x\  = '^TaijXj  i = 1,  2,  3.  (4.38) 

j= 1 

Since  the  transformation  must  preserve  the  distances  between  particles  in  the  rigid  body 
the  coefficients  form  a special  type  of  matrix  which  is  a representation  of  an  orthog- 
onal transformation.  This  matrix,  which  will  be  the  representation  of  an  operator  A,  is 
orthogonal  and  therefore  must  satisfy 


'y  ' Q’ijQ'ik  3jk-  (4-39) 

i 

The  operator  A may  be  thought  of  as  a change  of  basis  or  for  our  purposes  as  a trans- 
formation of  the  position  vector  r to  a different  vector  r': 

f = A r.  (4.40) 

It  is  important  to  note  that  a product  of  two  orthogonal  transformations  is  also  an  or- 
thogonal transformation.  If  we  restrict  ourselves  to  transformations  with  determinant 
+1  then  rotations  done  consecutively  are  equivalent  to  one  single  rotation. 
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4.3.2  Euler  Angles 

As  we  observed  in  the  previous  section,  rotations  can  be  mathematically  realized 
as  orthogonal  transformations.  Here  we  describe  how  these  transformations  are  con- 
structed and  practically  implemented.  The  simplest  rotation  is  one  about  a particular 
axis  through  and  angle  6.  For  such  a rotation  of  the  transformation  matrix  can  be  eas- 
ily calculated  using  basic  trigonometry,  for  example  a rotation  about  the  x-axis  in  three 
dimensions: 


/ 


\ 


V 


10  0 

A=  0 cos  6 sin  (9  • (4.41) 

0 — sin  6 cos  6 

In  order  to  describe  the  motion  of  rigid  bodies  in  the  mechanics  three  independent  pa- 
rameters are  needed.  The  choice  of  these  parameters  is  one  made  of  convenience  and 
appropriateness  for  particular  problems,  here  we  will  describe  two  of  the  more  popular 
formulations. 

The  most  popular  set  of  parameters  are  the  Euler  Angles.  They  are  defined  as  three 
successive  rotations  performed  in  a specific  sequence  through  three  angles.  Within  lim- 
its, the  choice  of  rotation  angles  is  arbitrary  but  the  most  common  is  to  rotate  the  system 
about  the  coordinate  axis  (as  in  Eqn.  4.41)  using  a particular  order.  One  such  sequence 
is  started  by  rotation  the  initial  system  of  axes,  xyz,  by  an  angle  (p  counterclockwise 
about  the  z-axis,  giving  the  new  axis  as  x'y'z.  In  the  second  stage  the  intermediate  axes 
are  rotated  about  the  x'-axis  counterclockwise  by  an  angle  9 to  produce  another  inter- 
mediate set  of  axes  x'y" z' . Finally  these  axes  are  rotated  again  about  the  x'-axis  in  a 
counterclockwise  direction  by  and  angle  ip.  The  three  angles  9,  (p  and  ip  constitute  the 
three  Euler  angles  and  they  completely  specify  the  orientation,  labeled  XYZ,  of  the 
new  system  relative  to  the  original  coordinates. 
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The  elements  of  the  complete  transformation  A can  be  obtained  by  writing  the  ma- 
trix as  the  triple  product  of  the  separate  rotations,  each  of  which  has  a relatively  simple 
matrix  form  similar  to  Eqn.  4.41 . The  complete  transformation 


A = BCD 


is  composed  of 


D 


^ cos  (f>  sin  (f)  0 ^ 

— sin  (f)  cos  (j)  0 

0 0 1 


and 


B 


1 0 0 
0 cos  9 sin  9 
^ 0 — sin  9 cos  9 

^ cos  'ip  sin  ■0  0 ^ 

- sin  ip  cos  ip  0 

0 0 1 


(4.42) 


(4.43) 


(4.44) 


(4.45) 


4.3.3  Cayley-Klein  Parameters 

Although  we  have  seen  that  only  three  independent  quantities  are  needed  to  specify 
the  orientation  of  a rigid  body,  there  are  occasions  when  it  is  desirable  to  use  more  than 
the  minimum  number  of  quantities.  The  Euler  angles  are  difficult  to  use  in  numerical 
computation  because  of  the  large  number  of  trigonometric  functions  involved,  and  the 
four-parameter  representations  are  much  better  adapted  for  use  on  computers.  We  can 
write  the  matrix  A in  terms  of  four  real  parameters  e0,  ei,  e2, 
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The  matrix  representation  is  then 


A = 


4 4-  e\  — e\  — e\  2 (eie2  + ^063)  2 (eie3  — eoe2) 

2 „2  1 „2  „2 


eo  ei 


2 (eie2  — eoes) 

^ 2 (eie3  + eoe2)  2 (e2e3  — eoei) 


e2  — e3  2(e2e3  + eoei) 


eo  - ei  ~ e2  + 4 


(4.46) 


The  parameters  obey  the  relation 


2 , 2 , 2 . 2 1 

e0  + el  + e2  + e3  — 1 


(4.47) 


and  are  derived  from  the  rotation  formula  of  Rodrigues[32],  This  formula  represents 
a single  finite  rotation  about  an  axis  to  transform  the  body  to  its  new  orientation.  It  is 
given  by 

r*  = f cos  <E>  + n (n  ■ r)  [1  — cos  $]  + (rx  n)  sin  <f>.  (4.48) 
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Here  r is  the  original  vector  which  is  rotated  through  an  angle  $ in  a plane  defined  by  a 
vector  normal  to  the  rotation  axis  n as  in  Fig.  4.2.  The  rotation  formula  can  be  rewritten 


Figure  4.2:  Vector  diagram  for  derivation  of  the  rotation  formula 


using  the  following  relations: 


<f> 

eo  = cos  — 

2 

$ 

e = n sin  — (4.49) 

2 


where  e is  the  vector  with  components  e0,  ex,  and  e3.  Asa  result  we  have 


r en 


:|)  + 2e  (e  • r)  + 2 (r  x e) 


(4.50) 


49 


which  can  be  shown  to  be  equivalent  to  the  orthogonal  transformation  given  by  Eqn. 
4.46. 

4.3.4  Rate  of  Change  of  a Vector 

The  dynamics  of  the  rotational  motion  of  a rigid  body  is  best  understood  by  exam- 
ining infinitesimal  rotations  first.  Consider  some  arbitrary  vector  r involved  in  a me- 
chanical problem,  such  as  the  position  vector  of  a point  in  the  body,  or  the  total  angular 
momentum.  Usually  such  a vector  will  vary  in  time  as  the  body  moves,  but  the  change 
will  often  depend  on  the  coordinate  system  to  which  the  observations  are  referred.  For 
example,  if  the  vector  happens  to  be  the  radius  vector  from  the  origin  of  the  body  set 
of  axes  to  a point  in  the  rigid  body  then,  clearly,  such  a vector  appears  constant  when 
measured  in  body  set  of  axes.  For  an  observer  fixed  in  the  space  set  of  axes  finds  that 
the  components  of  the  vector,  vary  in  time,  when  the  body  is  in  motion. 

The  change  d t in  a time  of  the  components  of  a general  vector  f as  seen  by  an 
observer  in  the  body  system  of  axes  will  differ  from  the  corresponding  change  as  seen 
by  an  observer  in  the  space  system.  We  can  derive  a relation  between  the  two  differential 
changes  in  f based  on  physical  arguments.  We  can  write  that  the  only  difference  between 
the  two  is  the  effect  of  rotation  of  the  body  axes: 

(d0  space  = W body  + (dF>  rotation  ' (4.51) 

Consider  now  a vector  fixed  in  the  rigid  body.  As  the  body  rotates  there  is  of  course 
no  change  in  the  components  of  this  vector  as  observed  from  the  body  frame.  The  only 
contribution  to  (d r)space  is  then  the  effect  of  the  rotation  of  the  body.  But  since  the 
vector  is  fixed  in  the  body  system  , it  rotates  with  it  counterclockwise,  and  the  change 
in  the  vector  as  observed  in  space  is  that  given  infinitesimal  version  of  Eqn.  4.46  which 
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can  be  written  as 

(d r\otatlon  = d(l  x f.  (4.52) 

The  time  rate  of  change  of  the  vector  fas  seen  by  the  two  observers  is  then 


space 


+ lJ  x r , 


(4.53) 


where  C u is  the  angular  velocity  of  the  body  defined  by  the  relation  udt  = df2.  The 
vector  Co  lies  along  the  axis  of  the  infinitesimal  rotation  occurring  between  t and  t + d t, 
a direction  known  as  the  instantaneous  axis  of  rotation.  The  magnitude  of  u measures 
the  instantaneous  rate  of  rotation  of  the  body. 

The  angular  velocity  can  be  expressed  in  terms  of  the  Euler  angles  and  their  deriva- 
tives in  both  the  space  and  body  frames.  In  the  body  frame  the  general  infinitesimal 
rotation  associated  with  u can  be  considered  as  consisting  of  three  successive  infinitesi- 
mal rotations  with  angular  velocities  = (/>,  u>g  = 9 and  = ip.  In  the  body  frame  in 
Cartesian  coordinates  these  become 


c o'x  = (p  sin  6 sin  ip  + 9 cos  ip 
uj'y  = (p  sin  9 cos  ip  — 9 sin  ip 

uj'z  = (p  cos  6 + ip.  (4.54) 

In  the  space  frame  we  have 

u)x  = 9 cos  <p  sin  ip  + ip  sin  9 sin  <p 
ujy  = 9 sin  (p  sin  ip  — ip  sin  9 cos  (p 

ujz  = ipcos9  + <p.  (4.55) 


4.3.5  Rotational  Inertia,  Angular  Momentum  and  the  Euler  Equations 

In  the  center  of  mass  frame  of  a rigid  body  the  total  angular  momentum  is 
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(4.56) 


Q 


Introducing  the  angular  momentum  it  is  possible  express  Eqn.  4.56  as 


L = I vecuj 


(4.57) 


where  I is  the  rotational  inertia  tensor.  In  the  CM  we  may  express  the  rotational  inertia 
tensor  as  a matrix  with  elements 


It  is  possible  to  diagonalize  I by  solving  the  appropriate  eigenvalue  problem  yielding  a 
transformation  matrix.  The  eigenvalues-/i,  /2,  and  /3- are  commonly  referred  to  as  the 
principal  axes  of  inertia.  We  can  transform  the  coordinate  description  of  the  orientation 
of  the  rigid  body  so  that  the  I is  diagonal.  This  simplifies  the  relevant  equations  greatly, 
the  angular  momentum  can  be  written  in  coordinate  form  as 


Ijk  ^ , ma  (ra  • ^aj^ak)  • 


(4.58) 


a 


(4.59) 


A set  of  equations  know  as  Euler’s  equations  can  be  derived  for  the  rotational  motion 
about  a fixed  point  using  a direct  Newtonian  approach.  We  consider  either  an  inertial 
frame  whose  origin  is  at  the  fixed  point  of  the  rigid  body,  or  a system  of  space  axes  with 
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origin  at  the  center  of  mass.  In  these  situations  we  have 


space 


(4.60) 


where  N is  the  external  torque  on  the  system.  This  equation  can  also  be  expressed  in 
terms  of  the  body  frame  derivitives  by 


space 


LU  X L. 

body 


(4.61) 


We  can  then  substitute  and  the  resulting  in  the  Euler  equations.  Expanding  the  equations 
gives 


I\U  1 ~ ^2^3  (I2  — h) 

= Nu 

I2UJ2  ~ <^3^1  {h  ~ h) 

= n2, 

I3U3  — LO1U2  (Ii  — I2) 

= n3. 

(4.62) 

For  torque  free  motion  we  have  = 0 so 


I\U>1  — U>2^3  {h  ~ h)  > 


I2UJ2  — ^3^1  {h  ~ h)  > 
-^3^3  = UJ1U2  ( I\  — I2)  ■ 


(4.63) 
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4.4  Quantum  Mechanical  Theory  of  Rotations 

4.4.1  The  Rotationally  Invariant  System 

Consider  a problem  where  V (p,  <f>,0)  = V ( p ).  Then  a change  of  basis  to  spherical 
coordinates  allows  us  to  write  the  eigenvalue  equation  for  H as 


- h 2 (&_  ld_  1 d2  \ 

2 p + pdp  + p2  dcf)2  ) 


+ V (p) 


(P,  <t>)  = ElpE  (p,  <f>)  . 


(4.64) 


Here  p denotes  the  mass  of  the  system.  The  angular  momentum  operator  in  three  di- 
mensions can  be  formulated  using  infinitesimal  rotations  which  gives  us  the  components 


Lx  — YPZ  — ZP  y 
L y — ZP  x — XP  2 

Lz  = XP,-YPx.  (4.65) 

These  components  obey  the  commutation  relations 

[Li,Lj]  = ih  Y.  £ijkLk  (4.66) 

k=x,y,z 

with  i,  j and  k running  over  the  coordinate  indices  x,  y and  z.  Additionally,  is  the 
Levi-Civita  tensor  also  known  as  the  antisymmetric  tensor  of  rank  3.  With  this  operator 
it  is  possible  to  rewrite  the  rotational  Hamiltionian  in  operator  notation 

H=5T+yW 


(4.67) 
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where  / is  the  rotational  inertia  and  L2  is  the  Euclidean  norm  of  the  angular  momentum 
operator.  Clearly  we  have 


[H,  Lj]  = 0 
H,  L2]  = 0. 


(4.68) 


Consequently,  L2  and  its  components  are  conserved.  Now  we  see  that  we  can  find  a 
basis  common  to  H,  L2,  and  one  of  the  components  of  L,  since  the  components  do  not 
commute  (see  Eqn.  4.66). 

4.4.2  Eigenvalue  Problem  of  L2  and  L, 

Following  the  example  of  the  harmonic  oscillator  we  will  now  examine  the  the  eigen- 
value problem  of  L2  and  Lz.  Let  us  assume  that  there  exists  a basis  | a/3 ) common  to 
both  of  these  operators: 


L2 1 a/3 ) = a\  a/3) 
hz\a(3)  = (3\a(3) . 


(4.69) 


We  now  define  raising  and  lowering  operators 


L±  — Lj  i ii-iy 


(4.70) 


which  satisfy 


[L2,  L±]  — 3zhL± 


(4.71) 


and  of  course 


L2,  L±]  = 0. 


(4.72) 
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The  final  results  of  the  eigenvalue  problem  are 


L2  | lm)  — l (l  + 1)  h2  | lm)  l = 0,  1,  2,  . . . 

L2  | lm)  = mh2  | lm)  nn  = — /,  l — 1,  l — 2,  . . . , l 


(4.73) 


for  integral  values  of  the  angular  momentum.  The  complete  solution  of  the  eigenvalue 
problem  actually  includes  half  integer  angular  momentum  but  these  are  not  of  interest  at 
this  time.  The  parallel  construction  between  harmonic  oscillators  and  angular  momen- 
tum begins  to  break  down  here  where  the  eigenvalue  is  bounded  from  above  and  below. 
Expanding  the  operators  in  the  coordinate  basis  allows  us  to  determine  the  eigenvectors 
for  the  rotationally  invariant  problem.  The  general  solution  includes  the  radial  part  while 
the  angular  part  of  the  function  is  given  by  the  spherical  harmonic  functions 


This  solution  can  then  be  used  for  particular  values  of  the  potential  to  arrive  a the  specific 


ipEim  ( pA,0)  = RElm  (r)  Y™  (0,0). 


(4.74) 


solution. 


CHAPTER  5 

COHERENT  STATES  AND  ROVIB  RATIONAL  DYNAMICS 


The  minimum  definition  of  a Coherent  State  (CS)  given  by  Klauder[33]  is  that  they  are 
a set  of  normalized  elements  | A ) of  a Hilbert  space  H.  They  share  the  following  two 
properties: 

(1)  continuity,  the  states  are  continuous  functions  of  a parameter  set  A 


lim|A)  = |A0),  (5.1) 

A— >Ao 

(2)  resolution  of  the  identity,  there  exists  a positive  measure  d/j,+  > 0 for  which 

1 = J dp+  | A)  ( A | . (5.2) 

Importantly,  the  set  of  CS  form  an  overcomplete  set  of  vectors.  Hence,  the  closed  linear 
span  of  Coherent  States  is  the  Hilbert  space,  H.  Coherent  States  have  proven  very  useful 
to  many  areas  of  the  physical  sciences,  and  occur  in  great  variety. 

5.1  Quasiclassical  Coherent  States 

A vital  property  of  some  CS  is  their  quasiclasscial  (QC)  dynamics.  A state,  | ip), 
is  said  to  be  quasiclassical  when  the  expectation  values  of  the  position,  momenta,  and 
energy  satisfy  the  classical  Hamiltonian  equations  of  motion.  Mathematically, 

Xqc  = (lp\x  \lp),  Pqc=  {lp\p\lp),  Hqc  = (lp\  H \lp)  , (5.3) 


satisfy 


_ dHqc  . = dHqc 

dxqc  Pqc  dpqc  ' 


(5.4) 
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The  average  position,  momenta,  and  energy  of  the  quasiclassical  state  evolve  in  time  as 
the  position  and  momentum  of  their  classical  analogs.  It  is  important  to  note  that  this  is 
not  the  same  demand  as  the  semi  classical  limit,  when  h 0,  is  invoked.  There  is  no 
guarantee  that  a quasiclasscial  state  even  exists  for  a given  Hamiltionian.  A method  to 
investigate  if  a Hamiltonian  allows  quasiclassical  states  is  to  apply  Ehrenfest’s  theorem 
and  show  the  resulting  equations  reduce  to  those  of  the  classical  case.  In  this  manner, 
we  can  show  that  the  canonical  CS  are  quasiclassical. 


The  Glauber  states  are  considered  the  canonical  Coherent  States  and  are  quasiclas- 
sical. These  states,  that  we  will  define,  | a)  are  associated  with  the  harmonic  oscillator 
Hamiltonian  Hvib  — huj(a+a  + |),  where  uj  is  the  angular  frequency.  The  harmonic 
oscillator  creation  operators  can  be  expressed  in  terms  of  the  self-adjoint  position  x and 
momentum  p operators 


where  m is  the  oscillator  mass.  Our  complex  parameter  a can  be  expressed  in  terms 
of  the  real  parameters  of  average  position  xa  = ( a \ x \ a ) and  average  momentum 

Pa  = ( a | p | a ) 


An  expansion  in  terms  of  harmonic  oscillator  energy  eigenstates,  {|  n) ; n = 0, 1, ..}, 
exists,  and  is  given  by 


5.2  Vibrational  Coherent  States 


(5.5) 


(5.6) 


(5-7) 
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Their  probability  distribution  over  the  harmonic  oscillator  states  is  given  by 


(5.8) 


The  Coherent  state  expectation  value  of  the  vibrational  energy  is 


Evib  =<  H >=  Hu 


(5.9) 


making  the  vibrational  existing  energy  a2  = ^djL.  So  the  distribution  of  the  vibrational 
levels  are 


exp 


V hco  ) n\ 


(5.10) 


Fig.  5.1  is  an  Example  for  a 2 = = 3. 


Figure  5.1:  Distribution  of  Coherent  States  over  harmonic  oscillator  energy  eigenstate 

«2  = = 3 


Notice  that  the  canonical  CS  are  associated  with  the  Hamiltonian  of  a vibrating  sys- 
tem. It  has  been  demonstrated  that  an  a posteriori  vibrational  analysis  in  terms  of  these 
harmonic  oscillator  CS  can  be  utilized  to  obtain  vibrationally  resolved  differential  cross 
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sections  when  the  vibrational  energy  levels  of  a product  wave  function  are  approxi- 
mately equidistant. 

A similar  derivation  can  be  made  for  CS  associated  with  the  Hamiltonian  of  a force 
free  rotor,  and  these  also  have  the  quasiclassical  property.  Additionally,  CS  have  been 
constructed  which  combine  rotational  and  vibrational  modes.  Here  we  have  a powerful 
tool  to  analyze  molecular  dynamics,  and  motivation  to  develop  a method  of  obtaining 
the  relevant  quantities  from  the  classical  evolution  of  a rovibrational  molecule.  It  is  our 
intention  to  generate  a general  method  of  analyzing  the  classical  vibration  and  rotations 
of  molecular  systems  in  order  to  evaluate  their  wave  functions  in  terms  of  quasiclassical 
CS. 


5.3  Proposed  Rotational  Coherent  States 

A new  rotational  CS  was  derived  by  Morales,  et  al.[34,  27]  in  the  spirit  of  the  Janssen 
CS  formulation.  Morales,  et  al.  expressed  these  CS  only  in  terms  of  the  integer  rota- 
tional states  so  that  it  is  quite  appropriate  to  describe  molecular  rotors  within  the  Quasi- 
Classical-Single-Determinental  END  theory.  Although  these  derived  CS  do  no  exhibit 
an  exact  quasi-classical  behavior  their  departure  from  the  classical  equations  can  be  eas- 
ily remedied.  Morales  was  able  to  show  that  the  exclusion  of  the  half-integer  rotational 
states  from  his  construction  leads  to  exact  quasi-classical  rotational  CS.  The  formulation 
of  half-integer  only  rotational  coherent  states  is  an  open  research  problem. 

Rotational  Hamiltonian  and  Related  Operators.  The  Hamiltonian  for  a free  ro- 
tating molecular  system  can  be  written  as[35] 


T 2 

TT 

Hr  _ 


lrot 


Ll 

T ~~ — b ~~~ 


2 It.  2 L 2 1. 


(5.11) 


where  /*  are  the  moments  of  inertia  for  the  principle  axis  and  L,  are  the  body  fixed 
components  of  the  angular  momentum.  For  molecular  systems  the  values  of  the  orbital 
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angular  momentum  are  restricted  to  integer  quantum  numbers  appropriate  for  bosons. 
The  space-fixed  components  of  the  angular  momentum  J(  (J2  = L2)  are  also  important 
for  our  discussion.  Taken  together  the  components  obey  the  following  commutation 
relations 

[Jj,  Jj]  i^ijk Jfc;  [Lj,  Lj]  zejjTjLfc  (5.12) 


and 


(5.13) 


where  ttjk  are  the  components  of  the  Levi-Civita  tensor.  Observe  the  asymmetric  com- 
mutation relationship  for  t components  of  L*.  From  these  relationships  we  can  construct 
a complete  set  of  rotational  eigenstates  |IMK)  such  that 


L2  |IMK)  = I (I  + 1)  |IMK) , 1 = 0,  1,  2,  ... 

Lz  |IMK)  = K |IMK) , K = 0,  ±1,  ±2,  . . . ± I (5.14) 

Jz  |IMK)  = M |IMK) , M = 0,  ±1,  ±2,  . . . ± I. 


These  rotational  eigenstates  can  be  represented  in  space  by  the  Wigner  D functions 


0,  A I |IMK> 


21  + 1 
. 87T2 


Dm  (<fi,  9,  V'j  ) 


(5.15) 


with  the  arguments  being  the  Euler  angles. 

Since  we  can  immediately  observe  that  the  the  Hamiltonian  Eqn.  5.11  commutes 
with  J i and  L;  the  Hamiltonian  eigenfunctions  must  satisfy 


HrofT^  — E°ot^  fM, 

J2^fM  = I(I  + l)^fM,  1 = 0,  1,  2,  ...,  (5.16) 

J^im  = 


M = 0,  ±1,  ±2,  ...±I, 
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where  the  superscript  a is  a third  quantum  number  to  label  a rotational  eigenstate  along 
with  I and  M.  Additionally,  it  can  be  shown  that 


[Hrot,  = i £ (L*L*  + L*Li)  ■ (5-17> 

ZIl 

i 

We  now  have  established  the  necessary  relationships  to  represent  the  rotational  Hamil- 
tonian eigenfunctions  \DfM  in  terms  of  the  angular  momentum  eigenstates,  or  in  mathe- 
matical terms 

*im  = £4Mq  |IMK>  (5.18) 

k 

where  the  coefficients  are  to  be  determined. 

Example:  Spherical  Rotor.  In  the  case  of  a spherical  rotor  the  moments  of  inertia 
are  equal,  ISr  — Ix  = Iy  = Iz,  and  the  eigenvalue  problem  reduces  to 


TT  _ _ 

nrot  ~ 2 ISR  ' 

*?M  = - 


J2 

2 ISR 


|IMK) 


rplMa  rp I 1(1+1) 

^ rot  rot  2Isr  ’ 


(5.19) 


5.3.1  Irreducible  Representation 

The  CS  under  investigation  are  not  group-related  but  are  built  up  in  analogy  to  the 
construction  used  for  group-related  CS  so  we  discuss  their  group  structure  for  complete- 
ness. The  irreducible  representation  of  the  rotational  eigenstates  |IMK)  is  the  semidirect 
product  of  the  50(3)  x 50(3)  with  an  Abelian  group.  The  generators  of  the  50(3)  Lie 
groups  are  the  Lt  and  the  J,,  respectively.  Each  have  an  associated  Lie  algebra  given  by 
their  commutation  relationships  which  satisfy 


[Lj , Lj]  — i^ijk  Lfc 


(5.20) 
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and 


[Jj,  Jj] 


(5.21) 


The  generators  of  the  Abelian  group  belong  to  a the  family  of  operators  Txu  (A  = 
0,  1,  §,  . . . ; ixv  = 0,  ±|,  ±1,  . . . , ±A)  which  satisfy 


and 


0, 


(5.22) 

(5.23) 


E(TA'TA'  = 4v,  <5-24) 

/3 

E(T^),T^  = <5'25> 


These  operators  are  a generalization  of  the  regular  spherical  harmonic  tensor  operators. 
If  A = | is  used,  as  in  Janssen[36],  which  gives  four  generators  for  the  Abelian  group. 
Unfortunately,  these  operators  connect  integral  (boson)  and  half-integral(fermion)  states, 
for  molecular  physics  we  should  use  A = 1 to  restrict  considerations  to  the  bosonic 
states. 

With  A = 1 this  leaves  nine  generators  for  the  Abelian  group.  In  order  to  estab- 
lish among  the  generators  of  the  elements  of  the  irreducible  representation  Morales,  et 
al.  generalized  the  properties  of  the  regular  two-subscript  spherical  tensor  operators 
relationships  with  both  Lj  and  J The  generalized  relationships  for  arbitrary  A are 


L2,Ty  = uT^,  (5.26) 

[jz>  = n T^, 


(5.27) 
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and 


L±,T£„]  = [A(A  + 1)-^=f1)]^T^tJ), 
J±,TjJ  = [A(A  + 1)  =F  1)]>  TfWA)„. 


(5.28) 


(5.29) 


The  action  of  the  generalized  T*v  operators  on  the  angular  momentum  eigenstates  is 
shown  to  alter  the  indices  and  multiply  the  states  |IMK)  by  a factor.  Moreover,  it  has 
been  shown  that  the  |IMK)  do  in  fact  span  the  irreducible  representations  of  the  product 
group.  It  is  important  to  note  the  action  of  one  particular  generator  of  the  Abelian  group 
on  one  particular  angular  momentum  eigenstates: 


5.3.2  Construction  of  Coherent  States 

The  Perelemov[37]  prescription  to  construct  coherent  states  associates  complex  pa- 
rameters with  the  generators  of  the  irreducible  elements  of  the  product  group.  In  this 
case  his  construction  would  require  three  complex  parameters  x , y,  and  z,  one  for  each 
quantum  number,  and  the  set  of  CS  will  be  generated  by 


Unfortunately,  this  set  of  CS  does  not  exhibit  the  quasi-classical  properties  for  a 
rotor.  Janssen[36]  altered  this  construction  to  generate  a set  of  quasi-classical  rotational 
CS.  Morales  modified  Janssen’s  construction  even  further  to  generate  CS  for  the  bosonic 
systems.  His  proposed  CS  are  given  by 


i - 1 - 1)  = ( 


21  + 
21  + 


(5.30) 


| xyz)  ~ e3:J+ezL+e2/T-1-1  |000) 


(5.31) 


-i-i|000),  (5.32) 
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where  the  operator  I is  defined  by 


I |IMK)  = I |IMK) 


(5.33) 


and  the  function  / (I)  is 


/(I) 


1 2I2  + 1 
2I2  — 1 


(5.34) 


The  operator  I admits  a representation  in  terms  of  the  Schwinger  boson  operators  for 
the  angular  momentum[36].  It  is  then  shown 


00  (vf(DTl,  nV 

exp  (yf  (I)  Tl,.,)  |000>  = £ |000) 

1=0, 1, ... 


I! 


°0  yl 

£ twU-i-i) 


(5.35) 

(5.36) 


and 


00  zn 


e'L+  |I  - I - I)  = £ =14  il  - I - I) 

^ — n il- 


n— 0 

21  zn 


- — r{21(21  - 1) . . . (2I  - [n  - l])}3n!5  |I  - I - 1 + n) 


71  = 0 


n\ 


= £ 


yl  + K 


tli  (I  + K)!t 


(21)! 

(I  K)! 


|I  — I — K) 


(5.37) 


where  n = I + K has  been  used  form  the  second  to  the  third  line.  Exchanging  M — > K 
and  L+  J+  an  analogous  expansion  for  ezJ+  |I  — I — I)  was  obtained.  All  of  these 
results  are  still  valid  for  integer  and  half-integer  angular  momentum  eigenstates.  The  set 
of  CS  can  now  be  expressed  as 


I xyz) 


exp 


yy*{  1 + xx*)2{\  + zz*)2 


x 
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y'  / (2I);2 l 2 x 

ti+m,.i7i+k 

V |IMK)  (5.38) 

I!  2 

These  states  are  normalized  so  that  ( xyz  \ \ xyz ) = 1 and  they  are  nonorthogonal.  By 
construction  the  rotational  CS,  | xyz),  satisfy  the  first  condition  of  the  definition  of  a CS 
(see  Eqn.  5.1).  The  second  condition  requires  a positive  measure  that  is  a resolution 
of  unity,  Morales,  et.  al  and  Janssen  were  only  able  to  construct  a weaker  version  of 
the  second  condition  in  which  the  measure  does  resolve  the  identity  but  is  not  positive 
everywhere.  The  measure  they  obtained 

d n±{x,y,z)  = -^{4[(1 +ZX*)  (1  + zz*)f  (yy*f 

7r° 

— 8 [(1  + xx*)  (1  + zz*)f  yy*  + l}dx  d y dz  (5.39) 


which  gave 


IMK)  = J d/r±  (x,  y,  z)  exp  -^yy*  (1  + xx*)2  (1  + zz*)2 


x 


X*l+My*l Z*l+K  X 


[(21)! 


I!(I  + M)!(I  — M)!(I  + K)!(I  — K)! 


xyz).  (5.40) 


The  goal  of  the  construction  of  the  rotational  CS  has  been  to  develop  quasi-classical 
states.  Morales,  et  al.  were  able  to  show  that  the  rotational  coherent  state  parameters 
can  be  related  to  the  relevant  physical  parameters  for  rotations.  Evaluating  the  operator 
averages  and  examining  their  dynamical  properties  they  were  able  derive[38]  relation- 
ships between  the  parameters  and  physical  quantities. 
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Connecting  Operator  Averages  and  Physical  Parameters.  Connecting  the  coher- 
ent state  parameters  with  physical  ones  was  accomplished  by  employing  the  following 
relationships 


x = e ia  cot 
z = Q+llf  cot 


where  0 < a,  p < 2n  and  0 < /3,  6 < n.  The  parameter  y is  expressed  as 


!,  = r^sin2(f) 


(5.41) 


where  0 < r < oo  and  0 < 6 < 2n.  The  CS  then  becomes 


xyz) 


| a (3 ip9  <f>r)  = | 

exp(-0  E { 


X 


X 


e-iMa  cosI+M 


e+iKip  cosI+K 


,1/2 


x— =|IMK). 
v/P 


CS) 

[(2I)!]2 

I!  (I  + M)!  (I  — M)!  (I  + K)!  (I  — K)! 


1 

2 


(5.42) 
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Calculating  the  operator  averages  reveals  the  significance  of  the  physical  parameters,  we 
observe 


( CS  | Lx  | CS  ) = r cos  <p  sin  6, 
( CS  | Ly  | CS  ) = r sin  p sin  9, 
(CS  |LZ  |CS)  = rcosO, 

< CS  1 1 1 CS  > = r 


( CS  | Jx  | CS ) = r cos  a sin  /3, 
( CS  | Jy  | CS ) = r sin  a sin  /3, 
< CS  I JZ|CS)  = r cos  P, 


(5.43) 


and 

( CS  | L2  | CS ) = ( CS  | J2 1 CS ) = r(r  + 2).  (5.44) 


It  is  then  possible  to  interpret  these  parameters.  If  follow  that  the  parameter  r is  the 
angular  momentum  modulus,  the  pairs  of  angles  ip,  and  9\  and  a,  and  (3  are  the  azimuthal 
and  the  polar  angles  of  ( CS  | L | CS  ) and  ( CS  | J | CS ) vectors  in  the  body-fixed  and 
laboratory  frame,  respectively.  The  angle  7 determines  the  relative  orientation  of  the 
body-fixed  frame  with  respect  to  the  space-fixed  one.  Finally,  the  probability  for  the  CS 
to  be  in  a rotational  state  1 1 M K ) is 


IMK 


X 


' Mif 

I!(I  + M)!(I-M)!(I  + K)!(I-K)! 

21+M  / P \ ■ 21-M  ( @ 
cos  ^ — sin  — 


cos 


21+K 


^ 2 
'0> 


x exp(— r) 


I! 


sm 


[(21)! 


21-K 


(I!  (I  + M)!  (I  - M)!  (I  + K)!  (I  - K)! 
xg(I  + K)(l-#-K)exp(— r)^, 


(5.45) 

■p(I+M)(l -p)(I~M)  (5.46) 


(5.47) 
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where  p = cos2  (f)  and  q = cos2  (f).  The  probability  is  the  product  of  binomial 
distributions  in  p and  q,  and  a Poisson  distribution  in  r.  We  observe 


00  I I 

£ £ £ Pimk  = 1-  (5.48) 

1=0, 1, ..  M=— I, ..  K=— I, .. 

Evolution  of  the  Coherent  State.  The  time  evolution  of  the  proposed  rotational  CS 
was  analyzed  using  Ehrenfest’s  theorem[38]  applied  to  the  operators  L j gives 


— < CS  | Li  | CS ) = i ( CS  | [Hrot,  Li]  | CS ) 

1 


( CS  | £ enkirr-  (ljl*  + L*Li)  I CS ) ■ (5-49> 

3 2Ak 


A quasi-classical  definition  for  the  angular  velocity  u was  introduced  where 


UJi 


; CS  | Li  | CS ) 
A,- 


(5.50) 


and  using  the  previous  averages  the  equations  of  motion  can  be  obtained 


ly 


(Iy 

(4 

(4 


4) 


4) 


(5.51) 


These  equations  are  very  similar  to  Euler’s  equations[32]  of  motion  for  a torque  free 
rigid  body.  They  demonstrate  that  the  derived  CS  are  almost  quasi-classical.  In  the  limit 
of  high  angular  momentum  (r  oo)  these  equations  would  match  the  classical  result. 


CHAPTER  6 

VIBRATIONAL  ANALYSIS 


Prony’s  method[39,  40]  (PM)  is  a technique  for  modeling  sampled  data  as  a linear  com- 
bination of  complex  exponentials  that  has  been  widely  used  in  engineering  and  the  phys- 
ical sciences[41,  42,  43,  44,  45].  It  is  highly  regarded  for  its  accuracy  and  stability.  It 
seeks  to  fit  a deterministic  exponential  model  to  the  data.  The  standard  PM  is  used  only 
to  analyze  one  degree  of  freedom,  but  we  will  require  an  analysis  of  many  more  degrees 
of  freedom.  We  will  build  an  analogous  method  to  analyze  our  simultaneously  rotating 
and  vibrating  system,  but  first  we  present  the  modem  version  of  Prony’s  technique. 


Assuming  one  has  the  N complex  data  points  x[i\  generated  from  p exponential 
terms,  the  Prony  method  will  estimate  the  terms  using  the  model: 


for  1 < n < N,  where  T is  the  sample  interval  in  seconds,  Ak  is  the  amplitude  of 
the  complex  exponential,  ak  is  the  damping  factor  in  I/seconds,  fk  is  the  sinusoidal 
frequency  in  Hz,  and  9k  is  the  sinusoidal  initial  phase  in  radians.  In  the  case  of  real 
data  samples,  the  complex  exponentials  must  occur  in  complex  conjugate  pairs  of  equal 
amplitude.  If  p is  even,  there  are  p/2  damped  cosines.  Ifpisodd,  then  there  are  (p— 1)/2 
damped  cosines  and  a single  purely  damped  exponential.  We  can  rewrite  our  p-exponent 
discrete-time  function  more  concisely  in  the  form 


6.1  Prony’s  Method 


p 


Ake[(ak+2nifk)(n-l)At+i9k] 


(6.1) 


p 


x[n]  = J2  hkzk  l. 


(6.2) 


k= 1 
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where  zk  are  the  time  dependent  factors  and  hk  are  the  time  independent  factors. 

Ideally,  one  would  like  to  minimize  the  error  x2  = J2iL\  (x[i]  - x[i])2  for  the  N data 
values,  with  respect  to  the  {hk}  and  {zk}  parameters  and  the  p exponents  simultane- 
ously. Evidently,  this  is  a difficult  nonlinear  problem,  even  if  the  value  of  p is  known. 
Iterative  algorithms,  such  as  steepest  descent  procedures  or  Newton’s  method,  have  been 
devised  to  minimize  this  error  with  respect  to  all  the  parameters.  Unfortunately,  these 
algorithms  have  proven  to  be  computationally  expensive.  The  Prony  method  embeds 
the  nonlinear  aspects  of  the  exponential  model  into  a polynomial  factoring,  for  which 
reasonably  fast  solution  algorithms  are  available.  Observe,  that  the  equation  may  be 
expressed  in  matrix  form  as 


4 

z° 

z2 

. z° 

p 

1 

1 
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H 

i— 1 
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zl 

z2 
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h2 
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z r1 
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7p- 1 

■ zp 

hp 

x\p\ 

If  a method  can  be  found  to  separately  determine  the  {zk}  elements,  then  Eqn.  6.3  rep- 
resents a set  of  linear  simultaneous  equations  that  can  be  solved  to  yield  {hk}.  Prony’s 
keen  insight  allowed  him  to  develop  just  such  a method. 

In  order  to  fit  a p exponential  model  using  PM,  at  least  2 p data  samples  are  required. 
An  exact  fit  to  the  p exponential  parameters  {hk}  and  {zk},  can  be  obtained  with  2 p data 
samples.  The  method  has  been  extended  to  optimize  the  2 p parameters  when  more  data 
is  available  using  least  squares  fitting. 

The  key  to  the  separation  is  to  recognize  that  the  Eqn.  6.3  is  a solution  to  a homoge- 
neous, linear,  constant-coefficient  difference  equation.  In  order  to  find  the  form  of  this 
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difference  equation,  first  define  the  polynomial  <f>(z)  that  has  the  {z\, ...,  Zk } as  its  roots, 

Hz)  = n(z  ~ Zk)-  ^6-4^ 

k- 1 

Expressed  as  a power  series,  this  becomes 

(f>{z)  = t a[m]zp-m,  (6.5) 

m— 0 

with  complex  coefficients  {a [m]}  such  that  a [0]  = 1.  Shifting  the  index  from  n to  n — m 
and  multiplying  by  the  parameter  {a[m]}  yields 

v 

a[m]x[n  — m]  = a[m]  Y hkZ^m~l,  (6.6) 

fc=i 

which  is  valid  for  p + 1 < n < 2p.  Summing  over  similar  products  produces 

Y a[m]x[n  - m]  = Y hk  Y a[m\'4~m~\  (6.7) 

m—0  k= 1 m= 0 

Making  the  substitution  = z^~p~l zf~m , we  observe 

Y a[m]x[n  - m]  = Y hk4P'  Y a\m\zk~m  = °>  (6.8) 

m=0  k= 1 m= 0 

since  the  right-hand  summation  is  the  previously  defined  polynomial,  (f>(zk)  = 0. 
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The  p equations  for  the  values  of  {a[n]}  in  Eqn.  6.8  may  be  expressed  as  the  p x p 
matrix  equation 


1 

H 

.*T3  , 

S-2 

1 

i — 1 

a[i] 

x[p  + 1] 

x\p  + 1]  x[p]  ...  x{2] 

a[2] 

= — 

x\p  + 2] 

x[2p  — 1]  x[2p  — 2]  ...  x[p] 

. . 

x[2p\ 

This  equation  demonstrates  that  with  2 p complex  data  samples,  it  is  possible  to  deter- 
mine the  {hk}  and  {zk}  parameters.  The  complex  coefficients  {ak},  which  are  functions 
only  of  the  time-dependent  {zk},  form  a linear  predictive  relationship  among  the  time 
samples. 

In  summary,  PM  consists  of  three  distinct  steps.  First,  the  coefficients  of  polynomial 
{a/J  are  obtained  by  solving  Eqn.  6.9.  Second,  the  roots  of  the  polynomial  {zk}  are 
calculated.  Finally,  the  matrix  equation  Eqn.  6.3  can  be  solved  for  the  parameters  {hk}. 
Then  the  original  exponential  parameters  can  be  obtained  by  the  following  relations 

ak  = In  \zk\  /At  fk  = tan~l[Im{zk}/Re{zk}]/2nAt  (6.10) 

Ak  = \hk\  9k  = tan~1[Im{hk}  / Re{hk}]  (6.11) 

6.2  Generalized  (GPM)  Prony’s  Method 

Given  a set  of  nuclear  trajectories  for  a general  molecule  in  motion,  a method  anal- 
ogous to  PM  to  estimate  the  vibrational  and  rotational  parameters  will  be  presented.  A 
computer  algorithm  is  implemented  to  utilize  this  method  in  order  to  extract  normal  co- 
ordinates and  rotational  motion  from  ENDyne  molecular  simulations.  In  the  regime  of 
low  energies  and  small  vibrations,  it  is  possible  to  separate  the  rotational  and  vibrational 
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motion.  The  following  expression  is  a model  for  the  center  of  mass  frame  nuclear  po- 
sitions of  an  N atom  molecule  in  the  semi-rigid  approximation  for  which  we  consider 
vibrational  and  rotational  motion  to  be  decoupled 


t- 1 


A + £ A, 


jcje 


[2niQj  (t—  1)  At+ifij  ] 


3= 1 


(6.12) 


With  v = 1..N,  and  j = l..p.  The  time  step  between  data  points  is  given  by  At,  p 
represents  the  number  of  vibrational  modes  present.  Here  ~&t,v  is  the  position  of  the  ith 
atom  at  time  index  t,  which  corresponds  to  time  At  (t  — 1);  Ot  is  the  rotation  matrix 
where  O0  = T;  is  the  equilibrium  position  of  the  ith  atom.  The  ttj  and  pj  are  the 
jth  normal  frequencies  and  phase  respectively.  Additionally,  7*^-  are  the  normal  mode 
amplitudes  for  the  jth  mode  and  ith  particle;  and  the  Cj  are  the  normal  mode  weights.  In 
the  case  of  a classical  molecule  p is  constrained  by: 


3N  - 6 N > 3 


p < < 


3N  — 5 N <3 


If  we  assume  that  the  time  steps  At  are  small  so  that  the  rotations  are  infinitesimal 
from  time  index  t to  (t  + 1),  a useful  result  from  classical  mechanics  is 


A,  = i^,+At  = ^xA,„+6, 


“ VJ  3 dt 

3=1 


, (6.13) 


where  ~uft  are  the  infinitesimal  angular  velocity.  Rearranging  gives  us  the  expression 


tt,u  = ^ t,u  - -tit  X tUu  = Ot±  2mnj)e^iW-VAt+^ 

3= 1 


(6.14) 
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The  vibrational  part  of  this  variable  is  obvious,  but  it  is  modulated  by  the  rotational 
motion.  This  will  interfere  with  the  separation  vibrational  and  rotational  motion.  So  we 
choose  as  our  fundamental  variable: 


A*  = = t T'vjCji (6.15) 

3= 1 

It  proves  to  be  more  convenient  to  transform  the  variable  > {xt,q}  from  a list 

of  N vectors  to  a list  of  3 N scalars.  So  the  variable  q = 1..3 N.  The  xt:Q  variable  is 
precisely  the  form  prescribed  in  the  PM,  which  allows  the  Prony  procedure  to  extract 
the  desired  information  for  this  model  of  several  degrees  of  freedom.  For  our  case  of 
real  valued  data  there  is  a modified  version  of  Prony’s  method.  The  modified  method 
replaces  the  linear  prediction  error  with  the  conjugate  symmetric  linear  smoothing  error, 
expressed  as 


X 


2 

s 


'y  ] [t7qM*£m+7i,g  T 7l,g 

71  = 1 


2' 


(6.16) 


The  fundamental  variable  of  the  model  is  quite  similar  to  that  of  PM,  except  that  we 
seek  to  determine  more  than  the  smoothing  coefficients  {gq[n]}.  We  also  seek  to  deter- 
mine the  elements  of  the  rotation  matrix,  and  if  we  assume  the  rotations  are  infinitesimal 
then  it  is  equivalent  to  determining  the  angular  velocity  { u?7}-  Instead  of  determining 
the  smoothing  coefficients  by  solving  a matrix  equation  similar  to  Eqn.  6.9,  we  seek  to 
minimize  the  smoothing  error  Eqn.  6.16  relative  to  the  smoothing  coefficients  {gq[n}} 
and  the  angular  velocity  { ~uft}  . 

Once  the  smoothing  parameters  are  determined  we  proceed  by  forming  the  polyno- 
mial 
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Mz)  = 9q[p\z2p  + + 9q[  1]^P+1  + 9q[  4zP  1 + ...  + g*[p]  = 0.  (6.17) 


The  roots  of  this  polynomial  are  of  exactly  the  same  nature  as  presented  in  the  de- 
scription of  Prony’s  method  allowing  us  to  calculate  the  frequencies  f 1j  according  to 
Eqn.  6.10.  Using  these  roots  to  solve  for  the  following  matrix  equation  for  the  hqj  — 


Tq,jcj  (27riOj) 


1 

is* 

I-4  o 

4 

..  z° 

p 

*1,9 

4 

4 

z1 

. . 

hq,  2 

= 

X2,q 

i 

1 

zr1 

7p- 1 

..  zp 

hq,p 

xp,q 

(6.18) 


Tq  j can  be  recognized  as  the  inverse  of  the  orthogonal  transformation  to  the  normal 
coordinates 


3AT 

X!  ^q  J ’ Tq,k  — 

4,k  j,  k = 1..3N 

(6.19) 

q=l 

3 TV 

3 N 

= )-£~K.r 

(6.20) 

5=1 

5=1 

Employing  these  relationships,  we  can  determine  Cj , Tqj. 

6.3  Implementation 

The  software  implementation  is  displayed  in  Fig.  6.1  using  the  Class-Responsibility- 
Collaborator  diagram  model  [46],  it  is  meant  to  give  an  overview  of  the  variables  and 
their  interaction  with  each  other.  Fundamental  to  utilizing  the  GPM  scheme  is  the  min- 
imization of  Eqn.  6.16.  The  minimization  method  we  have  currently  implemented  is 
the  Fletcher-Reeves-Polak-Ribiere  [47]  conjugate  gradient  method.  This  requires  cal- 
culation of  the  gradient  of  the  function  to  be  minimized.  Additionally,  it  requires  that 
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the  user  supplies  reasonable  initial  values  of  the  undetermined  parameters  {gq[n]}  and 

6.3.1  Angular  velocity  and  Rotation  Matrix 

An  initial  estimate  of  the  angular  velocity  is  given  by  solving  the  classical  angular 
momentum  equation.  The  angular  momentum  is  given  by 

tt  = x mj5t,v).  (6.21) 

U=1 


The  moment  of  inertia  tensor  is  given  by 


N 

a=l 


3 

~ Ra,jRa,k 

1 


(6.22) 


We  solve  the  angular  momentum  equation  Lt  — and  obtain  the  infinitesimal 

rotation  vector  T$t.  Determination  of  the  angular  velocity  allows  us  to  calculate  the 
rotation  matrix. 

The  rotation  matrix  can  be  written  as 


Ot  = n Sr. 


T— 1 


(6.23) 


Where  St  is  the  infinitesimal  rotation  matrix  with  the  property  StRt,u  — ~Rt+i,v  and  is 
represented  in  the  lab  frame  by 


1 

-ujt3At 

(jOt2At 

At 

1 

—uJnAt 

— cat2  Af 

iOt\At 

1 

(6.24) 


77 


Unfortunately,  the  resulting  matrix,  0T,  is  not  orthogonal.  Since  we  are  interested  in  the 
inverse  of  the  rotation  matrix  we  redefine  the  rotation  matrix  so  that  it  is  orthogonal,  to 
ensure  a well  defined  process. 

It  is  possible  to  represent  finite  and  infinitesimal  rotations  by  the  rotation  formula 
presented  by  Goldstein  [32] 

= Itt,  cob*  + it  (it  ■ ~Ht,v)  [1  — cos<f>]  + (lit,v  x it)  sin<E>,  (6.25) 

where  $ is  the  angle  of  rotation  about  the  normal  to  the  plane  of  rotation  it.  We  can 
reformulate  Eqn.  6.25  into  a more  useful  form  by  introducing  a scalar  ea  and  a vector 
~t  defined  as 

— it  sin  ^ (6.26) 

<f> 

e0  = cos  — (6.27) 

$ = ±|ut|Af.  (6.28) 

Where  the  sign  ± is  given  by  the  right  hand  rule.  With  further  manipulation,  Eqn.  6.25 
can  be  rewritten  as 

%+ir  = %,v(e2Q  -el -e22-  e\)  + 2 lt(lt  ■ %tV)  + x ^)e„.  (6.29) 


The  elements  of  the  rotation  matrix  can  be  calculated  by  expanding  Eqn.  6.25.  Hence, 
the  resulting  infinitesimal  rotation  matrix  is  redefined  as 


ei  + ef 


2 (e2ei  + eoea)  2 (e3ei  — eoe2) 


St  = 


2 (e2ei  - e0e3)  e20 


el  + e2 


2 (e3e2  + eoei) 


2 (e3ei  + eoe2)  2 (e3e2  — eoei) 


eo  - ei  - e2  + e2 


3 


(6.30) 
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Then  the  rotation  matrix,  Ot,  is  orthogonal.  The  inverse  of  this  matrix  is  a rotation 
through  an  angle  — <f>,  the  transpose: 


e0  — > eo 

(6.31) 

~i  — 

(6.32) 

el  + e\-e22-  e\ 

2 (e2ex  - e0e3) 

2 (e3ei  + eoe2) 

2 (e2ei  + eoe3) 

p2  _i_  ^2  ^2 

e-i  -t  e2  ^3 

2 (e3e2  — eoei) 

(6.33) 

2 (e3ei  — eoe2) 

2 (e3e2  + eoei) 

eo  - ei  - 4 + 4 ) 

expressing  at  in  terms  of  e0  and  "Ogives 


at  = 


2 arccos  (e0) 
- eg 


(6.34) 


6.3.2  Smoothing  Coefficients 

The  initial  estimate  of  the  smoothing  coefficients  is  given  in  matrix  form 


R-2plt  2p 


2 X2P 
Ox 


V 0p 


(6.35) 
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where  the  matrix  R2p  and  vector  l?2p  are 


r2p[0,  0] 


Rip  — 


r2p[2p,  0] 


r2p[0,  2p] 

r2p[2p,2p] 


(6.36) 


and 


^[1] 

1 

9*1 1] 


(6.37) 


\ 9*\p] 

The  elements  of  R2p  are 

M 

r2p[j,  k]  = ^ (xn_jq  ■ Xn-k}q  + Xn-p+j^q  • Xn_p+k,q)  . (6.38) 

n=2p+l 

A fast  algorithm  has  been  provided  by  Marple  [39]  which  can  solve  these  equations 
using  the  symmetric  covariant  method. 


6.3.3  Gradient 

For  the  minimization  of  Eqn.  6.16  we  have  chosen  the  Fletcher-Reeves-Polak- 
Ribierre  conjugate  gradient  method  [47]  which  requires  that  we  calculate  the  gradient 
with  respect  to  the  adjustable  parameters.  A few  intermediate  results  that  we  will  find 
useful  are 

or1  = It  s,-V+i. 

T—  1 


(6.39) 
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ddr1  '-1 


de 


T = 1 


= 0 


/ T=/  + l 


for  l > t 


(6.40) 

(6.41) 


The  functional  dependence  of  ~$tv  on  the  parameter  rit  is  given  by  the  following 
derivation 


— Ot~i  (V  + V t)  i 


(6.42) 


where  Vt  represents  the  vibrational  terms  of  the  model.  Substituting,  we  obtain  the 
following  result 

(6.43) 


ft  = 7?1+1. , 7?‘  - X 7t, 


Af 


O,  (1  + ? t+ 1 ) - ot_!  (V  + Vt) 


At 


— ~uft  x Ot- 1 


The  properties  of  the  infinitesimal  rotation  matrix  yields  the  following  result 


(6.44) 


Ot  {V  + Vt)  = Afrit  x Ot-!  {V  + Vt)  + Ot-!  (V  + Vt)  ■ (6.45) 


Employing  this  result  gives  us 


and 


ft  = 0 


dVt 


dt 


rit  = ()-xrit  = 


dVt 

dt  ' 


Consequently, 


rit  = rit  {riitri2, , ..rit,rit{rit)) 


(6.46) 


(6.47) 


(6.48) 
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So  the  gradient  of  the  prony  like  variable  is  given  by 


de^i  de _ 

where  = e^u^). 

The  first  set  of  derivatives  of  the  smoothing  error  are 


(6.49) 


9x 


dg„  [A:] 


= -2E 


M-p  ( p \ 

^ m,u  ^ m+n,u  T 9u  [^]  m— n,u\  J " T m—k,i> 

m=p+ 1 V 7i=l  / 

(6.50) 


for  A;  = l..p  and  v = 1..N.  The  second  set  of  derivatives  are 

9x1 


de 


Hi 


2L 


M— p / 

V \ 

— 'y  ] 9v[^\~^ m+n,v  + Si/[^]  ^ m—n,v  I 

m=p+l  \ 

71=1  / 

(6.51) 


/dV  p 

\ k= l 


m+k,v 


de^l 


+ s£[fc] 


cpfr 

v 'Xj  m— 


m—k ,v 


de^l 


(6.52) 


Prony  Analysis 
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Figure  6.1:  CRC  Diagram  for  Generalized  Prony’s  method  subroutine 


CHAPTER  7 

VIBRATIONAL  CALCULATION  RESULTS 


In  order  to  investigate  the  efficacy  of  our  GPM  we  applied  it  to  an  interesting  nontrivial 
system.  Using  ENDyne[8]  we  calculated  the  dynamical  evolution  of  water  for  various 
magnitudes  of  vibrational  excitation.  Our  study  of  water  was  motivated  chiefly  by  the 
large  amount  of  experimental  data  available  for  this  chemically  and  biologically  impor- 

Normal  Modes  of  Water 


tant  molecule.  It’s  three  atoms  and  bent  geometry  provide  a theoretical  model  which  is 
not  simple,  but  also  not  very  complex.  The  bent  shape  of  the  H20  molecule  is  a result 
of  the  character  of  its  molecular  orbitals[48]. 


Ligure7.1:  H20  Modes 
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Normal 

Freq. 

Sym. 

Mode 

±1  cm-1 

Species 

V\ 

3656.65 

ai 

1594.59 

ax 

3755.79 

bl 

Table  7.1:  Experimental  vibrational  frequencies  for  H20.  Source:  Shimanouchi[49] 


h2o 

V\ 

^2 

Expt 

3656.65  cm-1 

1594.59  cm-1 

3755.79  cm”1 

SCF 

+332 

+102 

+344 

CISD 

+ 135 

+44 

+ 147 

MBPT(2) 

+77 

+ 15 

+112 

SDQ-MBPT(4) 

+89 

+33 

+107 

CCSD 

+80 

+34 

+98 

CCSD(T) 

+73 

+26 

+92 

Table  7.2:  Values  for  the  vibrational  frequencies  using  different  levels  of  theory.  The- 
oretical numbers  indicate  deviation  from  experiment  ie  SCF=Expt  + Value.  Source: 
Bartlett[50] 


7.1  Vibrationally  Excited  H20 

The  3N  — 6 = 3 (3)  — 6 = 3 normal  modes  of  water  are  visually  represented  by  Fig. 
7.1  where  the  symmetric  stretching  mode,  v2  the  bending  mode,  and  u3  the  asym- 
metric stretching  mode.  The  current  experimental  values  are  given  in  Table  7.1  and 
the  theoretical  harmonic  vibrational  frequencies  given  by  various  electronic  structure 
calculations  for  these  modes  in  Table  7.2.  It  is  important  to  realize  that  both  theory  and 
experiment  do  not  represent  the  situation  we  are  presented  with  in  the  dynamics  of  nu- 
clei in  molecules.  These  calculations  rely  on  the  assumption  that  the  nuclei  experience 
a truly  harmonic  potential,  specifically  that  the  potential  they  encounter  is  proportional 
their  displacement  from  equilibrium.  In  real  molecules  this  can  only  be  considered  an 
approximation  especially  as  the  vibrations  lead  to  large  displacements  from  equilib- 
rium. The  effect  of  this  breakdown  of  the  harmonic  assumption  is  commonly  called 
Anharmonicity  and  is  visually  represented  in  Fig.  7.2. 
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Two  Potentials 


Figure  7.2:  Comparison  of  a truly  harmonic  potential  and  a representation  of  a realistic 
intramolecular  potential 

Several  important  features  of  the  realistic  intramolecular  potential  Vr  should  be 
noted.  Firstly,  near  the  equilibrium  distance  the  differences  with  the  harmonic  poten- 
tial are  small  breaking  down  at  greater  displacements.  Importantly,  unlike  the  harmonic 
model  Vr  exhibits  a disassociation  energy  De  at  which  the  intramolecular  bond  will 
break[5 1 ].  Also  the  potential  Vr  exhibits  an  asymmetry  about  the  equilibrium  distance. 
All  of  these  features  represent  a breakdown  of  the  harmonic  potential.  The  deficien- 
cies of  the  harmonic  approximation  posed  a problem  for  our  analysis.  This  means  that 
realistic  intramolecular  potentials  deviate  from  the  harmonic  potential  more  for  higher 
vibrational  excitations,  and  we  are  faced  with  the  fact  the  molecular  motion  will  not 
consist  strictly  of  harmonic  vibrations,  especially  for  highly  excited  molecules. 

Often  realistic  potentials  are  examined  by  using  empirical  and  semi-empirical  ap- 
proximations such  as  the  Lennard-Jones  potential  or  the  Morse  potential.  We  chose  the 
Morse  potential  to  represent  the  intermolecular  bond  stretching  potential  of  the  F[20 
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which  is  given  by: 


U = Deq 


(7.1) 


In  Fig.  7.3  we  see  an  abstract  representation  of  the  Morse  potential  and  its  associated 
quantum  mechanical  energy  levels.  Here  we  again  see  the  effects  of  the  realistic  poten- 


Morse  Curve 


Intemuclear  distance  (Req) 


Figure  7.3:  Morse  Potential:  Energy  Levels  for  the  Morse  potential 


tial  compared  to  the  harmonic  potential.  Since  the  potential  allows  for  bound  states  the 
energy  levels  are  discrete  but  the  energy  levels  are  not  evenly  spaced  as  they  are  for  the 
harmonic  approximation.  In  fact,  the  spacing  between  sequential  energy  levels  becomes 
increasingly  smaller  for  larger  energy  states.  Additionally,  the  disassociation  energy 
limits  the  number  of  energy  levels  as  opposed  the  infinite  number  of  levels  available  to 
for  the  harmonic  oscillator. 
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We  used  one  method  of  obtaining  an  expression  for  the  energy  levels  by  employing 
the  Dunham[52]  expansion  of  the  morse  potential.  The  Dunham  expansion  of  vibra- 
tional energy  levels  is 

°°  / 1 \l 
£„  = EOfn+j) 

= ("  + 5)  ('  " ("  + j)  x'j  IjJ,:  + " ' <7-2) 

where  uje  is  the  frequency  given  by  the  harmonic  potential  given  by  the  quadratic  term 
in  the  Taylor  expansion  of  Eqn.  7.1  and 


Xe  = 


Lde 
4 D~e 


(7.3) 


is  referred  to  as  the  anharmonic  constant.  Immediately  we  can  identify  the  frequency 
of  vibration  given  by  Morse  potential  by  making  an  analogy  with  the  structure  of  the 
harmonic  energy  levels.  The  Morse  frequency  is  given  to  first  order  by 


Xe  We- 


(7.4) 


Unlike  the  equilibrium  harmonic  frequency  coe,  the  Morse  frequency  decreases  with  in- 
creasing quantum  number  n and  the  difference  in  the  between  the  frequencies  of  neigh- 
boring energy  levels  also  decreases. 

Employing  conservation  of  energy  it  is  possible  to  derive  the  first  order  correction 
of  the  intermolecular  motion 

q = (Sin  + ^ C0S  (2cUet))  ' (7-5) 
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This  result  shows  that  the  first  order  correction  of  the  harmonic  motion  given  by  the 
Morse  potential  generates  motion  which  is  simply  a sum  of  sinusoidal  functions.  This 
is  precisely  the  type  of  signal  required  to  perform  the  Prony  analysis.  Since  we  are  only 
concerned  with  the  low  lying  vibrational  states  we  will  leave  the  discussion  of  higher 
order  corrections  for  a later  time.  But  they  are  also  expressible  as  sums  of  sinusoids. 
Examining  Eqn.  7.5  reveals  that  the  effect  of  the  Morse  potential  on  the  displacement 
of  the  intermolecular  distance  for  a particular  vibration.  The  second  term  | cos  (2 uet) 
oscillates  with  twice  the  fundamental  frequency  of  the  main  component  of  the  motion. 
As  a result  the  motion  will  appear  assymetric  with  the  peaks  of  the  main  component 
being  more  spread  out  and  the  troughs  being  more  narrow  as  in  Fig.  7.4. 

Qualitative  Comparison  of  Harmonic  and  Anharmonic  Motion 

Anharmonic  Motion  to  First  Order 


Harmonic 

Anharmonic 


Figure  7.4:  Comparison  of  the  motion  in  a harmonic  potential  versus  the  first  order 
correction  due  to  an  anharmonic  potential 
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Armed  with  this  knowledge  of  the  effect  of  anharmonicity  on  the  vibrational  motion 
of  a molecule  we  then  investigated  the  ability  of  GPM  to  analyze  realistic  molecular 
motion.  We  did  this  by  performing  a series  of  Electron  Nuclear  Dynamics  trajectories 
for  only  H20  with  a single  vibrational  mode  excited  to  varying  degrees.  The  results 
of  the  trajectories  of  the  symmetric  bending  mode  of  water  from  equilibrium  are  given 
in  Fig.  7.5.  Here  the  effects  described  previously  are  plainly  seen.  The  frequency 
of  vibration  decreases  with  increasing  excitation  in  energy  and  there  is  a noticeable 
difference  between  the  peaks  and  troughs  of  the  motion. 


Vibration  of  Water 


Figure  7.5:  Endyne  Trajectories:  Excitation  of  the  stretching  mode  of  water  w/  fre- 
quency 4140  cm-1  using  STO-3G 

In  order  to  account  for  the  anharmonicities  we  applied  the  GPM  with  twice  as  many 
anticipated  signals;  which  was  motivated  by  the  observation  of  the  appearance  of  addi- 
tional sinusoidal  terms  predicted  from  our  Morse  potential  model.  Applying  the  GPM 
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in  this  way  to  these  resulted  in  the  data  in  Table  7.3.  The  additional  sinusoidal  terms 
were  clearly  measured  by  our  GPM  and  found  to  have  magnitudes  as  predicted  by  Eqn. 
7.5.  Moreover,  the  frequencies  decline  with  increasing  energy  as  predicted  by  Eqn.  7.4. 


ttyib 

hui 

loo  (cm  x) 

Ui  (cm  x) 

ratio 

.05 

4137 

8274 

2.0000 

0.1 

4134 

8269 

2.0002 

0.5 

4115 

8233 

2.0007 

2.5 

4006 

8194 

2.0500 

5.0 

3920 

7833 

1.9987 

Harmonic  :cue 

4140 

— 

— 

Table  7.3:  GPM  Analysis:  Predominant  frequencies  calculated  by  GPM  for  ENDyne 
trajectories  of  varying  excitations 


We  used  this  decrease  to  calculate  the  disassociation  energy  through  the  anhar- 
monic  constant.  Performing  a linear  regression  using  the  data  from  Table  7.3  (see  Fig. 
7.6)  we  obtain  a value  of  xe  ~ -012  which  corresponds  to  a disassociation  energy  of 
De  « .41  Hartree  or  n « 20.  It  should  be  kept  in  mind  that  these  results  are  the  product 
of  a first  order  correction  and  that  the  disassociation  energy  corresponds  to  a physical 
parameter  far  away  from  the  equilibrium.  We  observed  that  these  values  demonstrate 
that  our  morse  potential  approximation  is  consistent  with  physical  reality  and  provides  a 
useful  tool  in  analyzing  the  effects  of  anharmonicity  on  our  GPM.  Importantly  the  anal- 
ysis demonstrates  that  the  anharmonicities  do  not  cripple  our  ability  to  use  the  methods 
we  developed. 

7.2  Prony’s  Method,  Condition  Numbers  and  Numerical  Stability 

Our  development  of  Generalized  Prony’s  Method  is  for  the  purposes  of  numerical 
calculation  and  analysis.  Whenever  computers  and  numerical  methods  are  used  we 
must  remain  conscious  of  numerical  instability  and  loss  of  precision.  In  the  case  of 
the  Generalized  Prony’s  method  we  are  presented  with  two  sources  of  computational 
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Fit  of  the  Anharmonic  Constant 


Figure  7.6:  Best  fit  for  the  anharmonic  constant 


interest.  The  first  is  the  solution  of  a polynomial  (j>  (z),  but  here  we  focus  on  the  second 
which  is  the  solution  of  the  matrix  equation  given  by 

(zHz)  x — ZHb.  (7.6) 


where  Z is  a matrix  with  Vandermonde  structure[39].  This  matrix  Z is  a function  of  the 
roots  of  cj)  (z) 


(7.7) 


where  p is  the  number  of  exponentials  in  the  Prony  signal  as  well  as  the  order  of  z in 
4>  and  n > 2p  is  the  number  of  data  points  used.  The  solutions  of  (j)  are  related  to  the 
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frequencies  present  in  the  signal  by 


Zk  = exp  2ni£lk  At , 


(7.8) 


where  f)  is  the  frequency  and  At  is  the  sampling  period.  Of  course  we  are  primarily 
interested  in  the  matrix 


( 

7n 


7ip 


\ 


A = ZHZ  = 


y Tpi  ' ' ' 7pp  y 


where 


7 jk  — Ikj 


= E («?*)" 

m= 0 


n 


if  ^ 1 

if  = 1 


(7.9) 


(7.10) 


7.2.1  Method 

The  numerical  stability  of  the  solution  of  Eqn.  7.6  is  related  to  the  condition  number[47]. 

In  its  simplest  terms  the  condition  number  is  a measure  of  the  sensitivity  of  the  solution 

— * — * 

of  a linear  algebraic  system  Ax  = b with  respect  to  changes  in  vector  b and  in  matrix  A. 

It  is  calculated  as  the  product  of  the  norm  of  A and  the  norm  of  A"1  or  More  explicitly 


K = 


A||||A 


(7.11) 


The  norm  ||A||  in  this  instance  is  given  by  the  maximum  of  the  sum  of  magnitudes  of 
each  row,  or  formally 


max 


N 


E IAj 


(7.12) 


1 <j<N 
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The  main  computational  importance  of  the  condition  number  is  as  an  estimate  of 
the  precision  required  to  solve  the  algebraic  system.  An  ill  conditioned  matrix  is  one 
whose  condition  number  is  greater  than  or  of  the  order  of  the  reciprocal  of  the  machine 
precision,  i.e.  106  for  single  precision  and  1012  for  double  precision.  We  will  consider  a 
matrix  to  be  poorly  conditioned  if  the  square  of  the  condition  number  is  greater  than  the 
reciprocal  of  the  machine  precision. 

While  a general  result  relating  the  condition  number  to  the  parameters  of  our  par- 
ticular implementation  is  desirable,  we  first  only  investigate  our  application  of  GPM  to 
the  H20  molecule.  This  molecule  has  three  normal  mode  frequencies.  For  the  Prony 
analysis  this  indicates  that  the  number  of  exponential  parameters  is  six,  consisting  of 
three  complex  conjugate  pairs.  The  free  parameters  in  this  analysis  are  the  sampling 
period,  At,  and  the  number  of  data  points  n. 

We  will  evaluate  the  condition  number  for  several  reasonable  values  of  these  param- 
eters. The  range  of  the  sampling  period  will  be  At  = [1  au,  320  au]  with  the  upper  limit 
chosen  to  be  half  of  the  period  of  the  smallest  angular  frequency  since  the  periodicity  of 
the  signal  makes  any  larger  values  redundant.  The  number  of  data  points  ranges  from 
n = [2 p,  50p]  where  p — 6 the  minimum  number  of  points  needed  to  complete  the  full 
Prony  analysis,  the  upper  bound  was  chosen  to  be  a reasonably  large  number  of  data 
points.  Additionally,  we  will  investigate  the  effect  of  small  errors  in  the  frequency  12. 

7.2.2  Results 

The  parameters  obtained  from  the  first  step  of  the  Prony  analysis  are  the  solutions, 
z,  to  the  polynomial,  </>.  The  location  of  the  roots  in  the  complex  plane  can  be  observed 
in  Fig.  7.7.  In  the  presence  of  noise  in  the  prony  signal  the  location  of  the  roots  no 
longer  lie  on  the  unit  circle  in  the  complex  plane[39].  For  simplicity  and  convenience 
of  our  argument  we  disregard  the  phase  of  the  roots  so  {zj}  expressed  in  terms  of  the 
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Distribution  of  Solutions  in  the  Complex  Plane 
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Figure  7.7:  Location  of  the  roots  of  </>  (z)  in  the  complex  plane 

angular  frequencies  are  given  by 

Zj  — exp  ifijAt, 

where 

= 4140  cm-1  Q2  = 2170  cm-1  fi3  = 4390  cm-1. 

The  conversion  factor  to  atomic  units  is  c = 219474  (cm  au)_1  resulting  in 
fii  = 0.01886  air1  Vt2  = 0.009887  an"1  = 0.0200  an-1. 


Along  the  Unit  Circle 


(7.13) 


(7.14) 


(7.15) 


95 


The  corresponding  periods  associated  with  each  angular  frequency  are  given  by 


Therefore  we  have 


(7.16) 


T\  = 333.2  an 


635.6  au 


T,  = 314.2  au 


(7.17) 


The  solutions  then  are 


z\  = exp  iQiAt  Z2  = exp  ifl2At  Z3  = exp  iQ^At 
Z4  — exp  — At  z5  = exp  —iil2At  = exp  — iQ3At 


(7.18) 


and 


Z = 


/yTl  /y  77  /yTl  /y  77/  /yTl  /yTl 

zl  z2  z3  z4  z5  z6  y 


(7.19) 


7.2.3  Discussion 

In  the  analysis  of  the  numerical  stability  we  varied  the  length  of  the  sampling  period 
versus  the  number  of  data  points  used  in  our  Generalized  Prony’s  Method.  The  measure 
of  that  stability  was  the  logarithm  of  the  condition  number  and  we  can  observe  the 
results  in  Fig.  7.8.  The  most  important  aspect  of  this  graph  is  that  the  great  majority 
of  the  area  indicates  that  the  matrix  A in  Eqn.  7.9  is  well  conditioned.  We  should 
note  that  this  graph  displays  a broad  range  of  values  for  the  independent  parameters 
approaching  the  extreme  values  for  both.  In  our  applications[53]  we  used  a sampling 
period  of  At  = 70  au  and  number  of  data  points  n = 15 p = 180  which  is  located  in  a 
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Logarithm  of  the  Condition  Number 
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Figure  7.8:  Color  map  displaying  the  logarithm  of  the  condition  number  as  a function 
of  the  number  of  data  points  versus  the  length  of  the  sampling  period 


well  conditioned  region  with  typical  condition  numbers  for  all  of  these  calculations  in 
our  recent  paper  ranging  from  n = 1 ~ 100. 

An  interesting  feature  of  Fig.  7.8  that  is  apparent  are  the  lines  of  high  condition  num- 
bers occurring  at  particular  values  of  A t for  example  those  near  A t & 160,  210,  320  an 
This  feature  can  be  easily  understood  by  calculating  the  exponent  of  z3  for  a period  of 
At  = 160  au  where  we  observe 


ift3At  ~ i 0.0200  x 160  ~ 7 xi.  (7.20) 

As  we  see  the  sampling  period  is  approximately  equal  to  the  half  the  period  of  the  signal 
component.  Consequently,  the  values  of  z%  « ( — 1)A  and  we  will  see  a similar  effects 
for  all  {zj}.  This  periodicity  will  create  a matrix  Z with  many  elements  which  are  nearly 
identical  in  many  rows  which  is  the  classic  case  for  a poorly  conditioned  matrix  which 
will  be  reflected  in  the  matrix  A.  Therefore,  we  observe  that  particular  values  of  At 


97 


will  lead  to  poorly  condition  matrices  due  to  resonance  between  the  sampling  period 
and  angular  frequencies  present  in  the  data. 

The  region  of  large  condition  number  lying  in  the  bottom  left  comer,  coincides  with 
small  sampling  periods  and  few  data  points.  It  occurs  for  a similar  reason  as  the  lines  we 
observed.  A small  sampling  period  used  for  the  case  of  such  large  angular  frequencies 
leads  to  the  situation  where  is  small.  Again  we  observe  that  the  matrix  Z will  have 
many  rows  that  are  approximately  equal  to  one  another,  particularly  if  the  number  of 
data  points  used  is  small.  We  also  observe  a greater  area  of  large  condition  numbers  in 
the  upper  right  comer  of  Fig.  7.8 

Finally  to  establish  the  validity  of  this  analysis  we  performed  the  same  calculations 
with  a 10%  adjustment  to  each  of  the  frequencies.  As  we  expected  the  qualitative  fea- 
tures were  the  same  and  a minimal  quantitative  effect  was  the  result. 

7.2.4  Conclusion 

We  can  conclude  from  our  analysis  that  presence  of  Vandermonde  matrices  in  Prony’s 
method  is  not  a limitation  for  our  analysis.  Yet  we  must  be  judicious  in  our  choice  of 
sampling  period  with  respect  to  the  angular  frequencies  present.  We  can  anticipate  that 
for  certain  combinations  of  frequencies  there  may  not  be  such  a large  region  of  well 
conditioned  linear  algebraic  systems.  Especially  in  the  case  of  very  small  angular  fre- 
quencies present  in  the  same  signal  as  larger  angular  frequencies.  In  those  cases  other 
techniques  may  have  to  be  applied. 

7.3  Application:  Vibrational  analysis  of  END  trajectories 

In  order  to  test  GPM  in  a realistic  setting, an  analysis  of  a trajectory  of  H+  + H20 
at  46  eV  with  an  impact  parameter  of  2.7  au,  generated  by  Hedstrom,  et  al.  using  the 
ENDyne  program  is  carried  out.  The  data  analyzed  consists  of  a time  series  of  the 
positions  of  the  various  atoms. 
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Figure  7.9:  Vibration  over  80  au  of  time  for  a low  energy  transfer  collision  of  H+  + H20 
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Figure  7.10:  Rotation  over  800  au  of  time  for  a low  energy  transfer  collision  of 
H+  + H20 

7.3.1  Prony  Analysis  of  H20 

Low  Energy  Transfer.  We  analyzed  the  motion  the  water  molecule  that  is  the  mov- 
ing dynamically  on  an  average  weighted  surface  of  H20  and  FI20+  after  the  collision 
took  place.  The  data  was  converted  to  the  center  of  mass  frame.  Observe  that  the  pro- 
jectiles minimum  distance  from  the  target  is  still  quite  large.  The  resulting  dynamics 
yields  a water  molecule  that  has  become  slightly  vibrationally  and  rotationally  excited 
which  can  be  observed  in  Fig.  7.9  and  Fig.  7. 10. 
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V\ 

^2 

^3 

freq  (cm  *) 

3562 

1595 

3756 

STO  - 3 G 

4140 

2170 

4390 

Prony 

4014 

2093 

4299 

E/yi  & (&u) 

.009 

.0015 

.0024 

^ vib 

hu 

0.44 

0.16 

0.13 

quanta(au) 

0.49 

0.25 

0.52 

Table  7.4:  Resulting  coherent  state  parameter  a calculated  rom  the  application  of  GPM 

We  applied  the  GPM  method  to  the  END  trajectory.  We  observed  that  the  Prony-like 
variable,  Xt}i}Cc,  exhibited  high  frequency  low  amplitude  noise.  In  order  to  reduce  the 
impact  of  the  noise  in  the  signal  two  important  techniques  were  used.  Firstly,  the  data 
used  for  the  prony  analysis  is  sampled  from  the  original  data  at  a sampling  rate  appro- 
priate for  the  lower  frequencies  of  interest.  Secondly,  the  number  of  sinusoidal  signals 
was  overestimated,  as  suggested  by  Marple[39],  to  be  6 modes  instead  of  3.  This  has 
the  effect  of  allowing  the  prony  procedure  to  accommodate  some  of  the  noise  without 
sacrificing  the  estimation  of  the  higher  amplitude  sinusoidal  signal  we  are  interested  in. 
An  example  of  the  data  analyzed  is  displayed  in  Fig.  7.11 

After  the  GPM  we  are  presented  with  a list  of  frequencies  and  amplitudes.  We 
judiciously  select  the  3 frequencies  with  the  highest  amplitude  to  distinguish  our  desired 
signal  from  the  noise.  The  resulting  frequencies  and  amplitudes  are  shown  in  Table  7.4. 
The  theoretical  modes  are  compared  to  those  calculated  numerically  from  the  GPM  in 
Fig.  7.12. 

We  used  the  frequencies  and  amplitudes  obtained  from  our  data  analysis  to  generate 
a quasiclassical  coherent  state  (QCCS)  corresponding  to  our  classical  trajectory.  The 
resulting  distribution  over  harmonic  oscillator  energy  eigenstates  is  displayed  in  Fig. 
7.13.  As  we  expected  the  molecule  is  not  very  vibrationally  excited  due  to  the  large 
impact  parameter  of  the  projectile.  Also,  the  frequencies  determined  from  GPM  are 
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Vibrational  Analysis  of  Water 


Figure  7.11:  Comparison  of  the  actually  intermolecular  distance  over  time  to  that  calcu- 
lated by  GPM 


Normal  Modes  of  Water 
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Scale  Drawing  of  Modes 
Obtained  from  Prony  Analysis 


Figure  7.12:  Comparison  of  the  normal  modes  of  water  compared  to  those  calculated  by 
GPM 


slightly  lower  than  those  obtained  from  the  basis  calculation.  This  is  likely  due  to  the  fact 
that  in  the  ENDyne  implementation  the  molecule  is  represented  on  an  average  surface 
between  water  and  the  ionized  H20+.  Also,  it  is  possible  that  it  is  the  result  of  the 
presence  of  anharmonicities  either  of  which  would  reduce  the  frequencies. 

High  Energy  Transfer:  Violent  H20  Trajectory.  In  Fig.  7.14  there  are  numer- 
ous frequencies  present  that  cannot  be  accounted  for  by  just  normal  modes.  At  first 
this  seems  to  indicate  that  GPM  is  exhibiting  spurious  numerical  answers.  Observation 
shows  that  the  extra  frequencies  occur  at  regular  intervals,  it  maybe  hypothesized  that 
the  extra  signals  may  be  due  to  anharmonicities  in  the  molecular  bonds  due  to  a break- 
down of  the  Semi-Rigid  Approximation.  Central  to  the  SRA  is  the  description  of  the 
intra-molecular  bonds  as  ideal  springs  or  more  precisely  as  quadratic  potentials.  In  real 
molecular  bonds  this  approximation  is  only  valid  for  extremely  small  oscillations  about 
the  minimum  of  the  potential. 
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Coherent  State  Distribution  for  Water 


Figure  7.13:  Distribution  of  the  coherent  states  over  the  the  harmonic  oscillator  for  the 
three  normal  modes  of  water  analyzed  by  GPM 
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atomic  distance  for  violent  h2o  trajectory 


Figure  7.14:  Intermolecular  bond  distance  for  violently  excited  water  molecule 

7.4  Anharmonicity 

The  analysis  of  ENDyne  trajectories  of  water  for  different  energy  transfers  demon- 
strated that  the  GPM  analysis  reports  more  frequencies  than  can  be  accounted  for  by  the 
Semi-Rigid  Approximation.  Real  bonds  exhibit  many  features  which  harmonic  poten- 
tials simply  cannot  take  into  account.  These  features  can  be  observed  in  the  following 
figure.  Observe  potential  energy  curves  of  real  bonds  are  asymmetric.  Also,  that  har- 
monic potentials  do  not  reflect  that  there  exist  a disassociation  energy  above  which  the 
bond  no  longer  exists.  As  the  energy  of  vibration  increases  the  harmonic  potential  be- 
comes a poor  approximation  to  the  real  potential. 

In  order  to  explore  the  effect  of  the  breakdown  of  the  SRA,  several  ENDyne  trajec- 
tories of  varying  excitation  energies  of  the  stretching  mode  of  water  were  performed. 
For  the  STO-3G  basis  this  particularly  mode  is  predicated  by  GAMESS[54]  to  oscillate 
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with  frequency  4140  cm-1.  In  the  following  figure  it  is  evident  that  as  the  energy  of 
oscillation  increase  two  important  changes  occur  in  the  character  of  the  motion. 

Firstly,  as  the  energy  increases  the  wavelength  of  the  oscillation  increases  which 
corresponds  to  a linearly  proportional  decrease  in  frequency.  Most  importantly  higher 
energies  the  molecular  vibrations  begin  to  lose  their  sinusoidal  character  due  to  the 
asymmetry  of  the  potential.  During  the  contraction  of  the  bond  the  particles  move  more 
quickly  then  when  the  bond  is  elongated.  The  trajectory  exhibits  this  in  that  the  width 
of  the  troughs  become  increasingly  smaller  and  the  widths  of  peak  become  increasingly 
larger  as  the  energy  of  the  oscillation  increases. 

The  GPM  analysis  is  predicated  on  the  assumption  that  the  motion  of  the  molecules 
will  be  sinusoidal  in  character.  The  effect  of  the  asymmetry  of  the  potential  seems 
to  undermine  this  assumption.  Fortunately,  it  is  possible  to  derive  an  approximation 
to  realistic  potentials  which  still  allows  the  motion  of  real  nuclei  in  molecules  to  be 
described  by  sinusoidal  functions.  One  such  approximation  which  reflects  feature  of 
real  bonds  is  the  Morse  Potential 


The  solution  of  the  Schrodinger  equation  given  this  potential  yields  the  following 
equation  for  the  energy  eigenvalues: 


U = Deq  [l  -e-a^~r) 


(7.21) 


(7.22) 


(7.23) 


The  coefficients  C/  decrease  rapidly  providing  corrections  to  the  harmonic  oscillator 
energy  eigenstates.  The  first  coefficient  is  Cj  = u>e  and  the  next  and  most  significant 
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correction  C2  = Xe^e  where  the  Anharmonic  Constant  Xe  is  given  by 


Xe  = 


(jJe 

4 K' 


(7.24) 


An  immediate  consequence  of  this  expansion  is  the  following  result  for  the  frequen- 
cies associated  with  the  energy  eigenstates: 


^ n 


CJe 


(7.25) 


Observe  that  the  frequency  reduces  with  increased  n.  Using  perturbation  theory,  solving 
for  the  classical  motion  of  a particle  in  this  potential  yields  the  following  result  to  first 
order: 

q = (sin  (( unt ))  + ^ cos  (2a ;„£))  (7.26) 

uj2e  \ 8 ) 

Consequently,  the  particle  motion  can  be  approximated  by  a linear  combination  of  sinu- 
soids allowing  GPM  to  be  applied  to  these  motions. 

Applying  GPM  to  the  set  of  trajectories  of  excitations  of  the  stretching  mode  of  the 
water  the  analysis  detected  two  predominant  frequencies.  The  following  table  shows 
the  frequencies  detected  for  increasing  energies  and  calculates  the  ratio  between  the  two 
detected  frequencies.  These  ratios  are  in  excellent  agreement  with  the  first  order  motion 
given  by  Eqn.  7.26.  We  observed  that  the  fundemental  frequency  and  the  anharmonic 
correction  had  exactly  the  relationship  predicted  by  the  first  order  motion  as  in  Eqn. 
7.26. 

Plotting  the  fundamental  frequencies  reveals  that  they  decrease  linearly  which  was 
predicted  by  the  Dunham  expansion.  Performing  a least  squares  fit  gives  a value  for  the 
Anharmonic  Constant.  The  fit  estimates  that  Xe  ~ 012,  which  is  consistent  in  magni- 
tude given  by  Shimanouchi[49],  Additionally,  this  value  corresponds  to  a disassociation 
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^ vib 

hu 

i00  (cm  1 ) 

co i (cm  *) 

ratio 

.05 

4137 

8274 

2.0000 

0.1 

4134 

8269 

2.0002 

0.5 

4115 

8233 

2.0007 

2.5 

4006 

8194 

2.0500 

5.0 

3920 

7833 

1.9987 

Lde 

4140 

— 

— 

Table  7.5:  Calculation  of  anharmonic  effects  for  water  vibrationally  excited  at  various 
energies 


energy  of  De  « .41  Hartree  or  n « 20  which  is  an  overestimate  but  remarkably  good 
for  a first  order  correction. 


7.5  State  Resolved  scattering:  H+  + H20 

Classical  nuclear  trajectories  were  calculated  using  the  ENDyne  code  with  a mini- 
mal STO-3G  basis.  With  classical  nuclei  the  initial  orientation  of  the  target  needs  spec- 
ification. Using  the  C2i,  symmetry  of  the  water  molecule  a coarse  rotational  grid  of 
74  points  is  generated  from  30  target  orientations.  Approximately  thirty  trajectories  of 
varying  impact  parameters  are  calculated  for  each  target  orientation.  Each  trajectory  is 
evolved  to  2000  au  of  time  to  provide  a large  set  of  data  for  force  free  nuclear  trajectories 
for  vibrational  analysis.  GPM  is  then  applied  to  each  trajectory  and  the  corresponding 
QCCS  is  determined.  Then  for  each  orientation  of  the  target  the  deflection  function  is 
obtained  and  the  state  resolved  cross  section  is  determined.  Finally,  all  the  orientations 
are  rotationally  averaged  to  compute  the  final  state  resolved  cross  section. 

We  use  GPM  to  identify  the  three  normal  modes  for  each  trajectory.  The  nuclear 
trajectories  of  the  post  reaction  vibrationally  excited  water  molecule  are  extracted  from 
the  ENDyne  calculation.  We  denote  them  by  R „ [f]  where  v corresponds  to  the  O,  Hi, 
or  H2  and  t is  an  index  of  the  time.  Our  calculations  were  performed  with  a time  step  of 
A t = 1 au  which  is  sufficiently  small  to  consider  the  time  between  fj  — > ti+i  to  make 
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the  approximation 


D„  [f]  = jtK  [t] 


Ru  [t  + 1]  — Ru  [£] 
At 


(7.27) 


We  generate  an  estimate  of  the  infinitesimal  rotation  vector  ut  from  equation  for  angular 
momentum  leading  to  the  calculation  of  the  rotation  matrix  O and  its  inverse.  This 
allows  the  calculation  of  the  variable 


xv[t\  = O *[£]  Dv[t]  - Co  x Rv[t } 


(7.28) 


we  then  form  the  matrix  equation  given  by  Eqn.  6.35  to  obtain  an  estimate  of  the  smooth- 
ing parameter  gk.  Then  \2  is  minimized  for  the  parameters  gk  and  Cot.  The  smoothing 
parameters  are  used  to  define  the  polynomial  cf)  ( z ) which  is  then  solved  for  the  parame- 
ters zk- 

The  frequencies  {flj}  and  the  phase  {ipj}  are  obtained  by  employing  the  relation- 
ships given  in  Eqn.  6.10.  The  parameters  hqj  = TqjCj  (2n i£lj)  are  obtained  from  the 
matrix  equation 


1 1 

Zl  z2 


7p- 1 

zv 


1 

51 

1 

t — 1 

1 

hq,2 

= 

x[2] 

1 

73 

x\p\ 

(7.29) 


The  parameters  Tq  j are  recognized  as  the  elements  of  the  inverse  transformation  from 
Cartesian  to  the  normal  coordinates  and  therefore  satisfy 


3N 

,jTq,k  = T/fc. 

9=1 


(7.30) 


Consequently,  employing  Eqn.  7.30  the  following  result, 
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(7.31) 


determines  the  parameter  {cj}  , which  are  the  relative  weights  of  the  normal  coordinates. 
Once  the  {cj}  have  been  obtained  the  {Tqj}  . can  be  calculated  readily  from  hQj  — 
TqjCj  ( 2iriQj ) . As  stated  the  { Tq can  be  identified  as  the  elements  of  the  inverse 
transformation  from  Cartesian  to  normal  coordinates  so  the  GPM  can  be  used  to  obtain 
normal  coordinates  and  normal  frequencies  from  classical  nuclear  trajectories. 

In  the  case  of  our  water  trajectories  we  identify  the  three  modes  corresponding  to 
bending  vx,  symmetric  stretch  v2  and  asymmetric  stretch  u3.  Using  the  corresponding 
amplitudes  and  frequencies  we  calculate  the  vibrational  energy  of  each  mode  and  the 

o 

corresponding  parameter  a~  = where  j = 1,  2,  or  3 is  the  mode.  These  are  used 
to  generate  the  associated  quasiclassical  coherent  states  corresponding  to  that  trajectory. 
Since  we  can  consider  each  mode  to  be  an  independent  oscillator  we  have  three  sets  of 
coherent  states  of  the  form 


So  from  those  coherent  states  the  probability  to  be  in  a particular  vibrational  state  is 
calculated.  The  probability  for  charge  transfer  Pc  is  also  calculated  from  the  Mulliken 
population  of  the  outgoing  proton  in  the  output  of  the  ENDyne  code. 


(7.32) 


The  probability  to  be  in  an  harmonic  oscillator  eigenstate  n is  given  by 


(7.33) 


Once  the  probabilities  and  deflection  functions  have  been  calculated  the  differential 
cross  section  (DCS)  is  calculated  for  each  orientations.  By  including  the  probability 
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Deflection  Function 


to  be  in  a particular  vibrational  state  then  the  state  resolved  differential  cross  section 
(SRDCS)  may  be  calculated.  The  classical  state  resolved  DCS  is  given  by 


da  ^ Pb 

dQ  _M^)N,sin(0)lf 


(7.34) 


The  probability  to  be  in  a particular  vibrational  state  in  the  inelastic  collision  is  given  by 

p = pi,plpL  (i  - pcl 

The  DCS  and  SRDCS  are  both  calculated  using  the  semiclassical  uniform  Airy 
approximation[55]  in  order  to  remove  the  singularities.  The  SRDCS  for  all  orienta- 
tions are  then  rotationally  averaged  over  the  initial  target  orientations  in  order  to  obtain 
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SRDCS  for  comparison  with  experiments  where  there  is  no  preferred  orientation  of  the 
target  molecule.  The  average  is  achieved  by  performing  an  integral  over  the  Euler  angles 
given  in  Jacquemin,  et  al.[56] 

7.5.1  Experimental  DCS 

We  compare  our  results  to  the  experiments  performed  by  Toennies,  et  al.[57].  Their 
experiment  consisted  of  an  ion  beam  of  protons  crossed  with  a target  beam  of  H20.  The 
experiment  was  designed  to  analyze  two  possible  processes  either  inelastic  scattering  or 
charge  transfer, 


H+  + H20  - 

-4  H+  + H20 

(7.35) 

H+  + H20  - 

->  h + h2o+. 

(7.36) 

We  only  considered  the  inelastic  scattering  process  in  our  investigation.  For  both  pro- 
cesses the  scattered  protons  or  H atoms  are  detected  in  the  perpendicular  plane.  To  mea- 
sure the  H-atom  signal,  the  protons  are  completely  suppressed  by  an  electrical  repelling 
field.  The  count  rates  are  determined  either  as  a function  for  the  laboratory  scattering 
angle,  or  as  a function  of  the  flight  time  providing  information  on  the  state  selected  cross 
sections  at  a given  angle.  Relative  translational  energy  distributions  were  obtained  from 
the  TOF  spectra  by  transforming  the  time  scale  into  the  energy  scale.  This  energy  was 
plotted  against  the  net  change  of  the  relative  translational  energy  in  the  collision.  The 
experiment  was  performed  at  two  collision  energies  of  Ecm  = 27.0  eV  and  46.0  eV, 
TOF  spectra  and  angular  distributions  of  either  the  protons  or  H atoms  were  measured 
at  scattering  angles  scattering  ranging  from  0°  to  15°  and  0°  to  8°,  respectively.  Energy 
loss  distributions,  P (A E),  as  obtained  from  the  TOF  spectra  of  protons  released  in  the 
inelastic  process.  The  spectra  exhibit  significant  vibrational  peak  structures  which  are 
only  partly  resolved  with  respect  to  all  three  H20  fundamental  modes.  Toennies,  et.al. 
employed  a deconvolution  procedure  in  order  to  estimate  the  state-specific  transition 
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ENDyne  Nuclear  Trajectory  H++H20 

Classical  Nuclei  (46  eV) 


Figure  7.15:  Typical  Electron  Nuclear  Dynamics  trajectory  for  H+  + H20 

probabilities.  This  is  based  on  the  assumption  that  the  spectra  are  composed  of  indi- 
vidual vibrational  transitions,  starting  form  the  ground  vibrational  state,  which  all  can 
be  approximated  by  Gaussian  distributions  of  different  heights  but  of  equal  half  widths. 
The  rotational  excitation  was  estimated  to  be  50  — 100  meV  by  comparing  the  energy 
loss  spectra  to  atomic  scattering  of  Ne.  The  resolution  was  determined  to  be  compa- 
rable to  the  H20  bending  mode  energy,  hv2  — 198  meV.  Consequently,  this  mode  is 
only  observed  as  a shoulder  of  the  elastic  peak  in  all  the  experimental  spectra.  However, 
distinct  peaks  at  an  energy  loss  of  about  450  meV  appear  which  is  assigned  to  excitation 
of  the  symmetric  (v{)  or  asymmetric  stretch  (z43).  Since  the  energy  difference  between 
the  two  stretching  modes  is  only  13  meV  the  experiment  was  unable  to  resolve  these 
two  modes.  The  convention  was  made  to  denote  both  stretching  modes  by  v\  to  include 
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Figure  7.16:  Energy  loss  of  scattered  proton  in  Toennies,  et  al. [57] 

either  the  symmetric  or  asymmetric  stretch  vibrations  or  both  of  them. 

At  energy  losses  beyond  the  energy  of  the  v\  fundamental,  the  spectra  are  composed 
of  various  overtone  and  combination  transitions.  These  overtones  were  identified  in  the 
Gaussian  fit  curves  by  assuming  excitation  of  the  [z/l5  0,  0]  and  [vx,  u2)  0]  progressions. 
The  experimenters  identified  the  relative  heights  of  the  fitted  Gaussian  distributions  with 
the  transition  probabilities  from  the  initial  ground  vibrational  state  to  the  respective  final 
states.  These  probabilities  served  as  as  a basis  for  deriving  state-selected  quantities 
characteristic  of  the  H+  + H20  inelastic  collision  dynamics. 

The  branching  ratio  of  [zq,  0,  0]  and  [v\,  u2,  0]  progressions  were  determined  at  all 
available  angles.  They  showed  a gradually  increasing  importance  of  the  stretching  mode 
versus  the  bending  with  increasing  scattering  angle.  The  average  vibrational  energy 
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transfers  were  calculated  from  the  transition  probabilities.  From  the  measured  angu- 
lar distributions  of  the  scattered  protons  relative  state  selected  DCS’s  for  various  final 
vibrational  states  were  derived  using  the  relations 


das(i)  n^d<7(®) 

= Ps{l)~zr 


(7.37) 


(7.38) 


The  results  were  digitized  from  Fig.  7.16.  in  paper  by  Toennies,  et  al.[57]. 


7.5.2  Results  and  Discussion 

Our  aim  was  to  develop  numerical  techniques  to  obtain  vibrational  and  rotational 
spectral  information  from  a set  of  semiclassical  molecular  trajectories.  We  were  moti- 
vated by  our  desire  to  compare  our  simulations  using  the  ENDyne  program  with  exper- 
imental rovibrationally  resolved  DCS  of  molecular  scattering  dynamics.  In  this  paper 
we  have  presented  our  technique  to  employ  GPM  in  combination  quasiclassical  coherent 
states  to  obtain  state  resolved  cross  sections  from  classical  molecular  trajectories.  Here 
we  will  present  the  results  of  GPM  application  to  a single  trajectory  and  also  to  more 
stochastic  processes  such  as  the  calculation  of  state  resolved  DCS.  Additionally,  we  will 
compare  these  theoretical  calculations  with  those  of  Toennies  et  al.[57]  and  discuss  the 
dynamics  of  the  H+H20  inelastic  scattering  process. 

As  an  example  we  present  the  results  for  the  GPM  application  to  one  typical  H+  + H20 
trajectory  in  Table  7.6.  We  compare  the  ab  initio  results  for  the  three  normal  modes  with 
those  obtained  from  the  GAMESS  program[54,  58]  using  and  the  GPM  results.  First  we 
observe  that  the  the  GAMESS  results  were  obtained  using  the  STO-3G  basis  that  was 
used  in  the  ENDyne  calculations.  Consistent  with  common  experience  the  frequencies 
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generated  from  the  SCF  calculation  of  GAMESS  overestimates  those  frequencies  by 
approximately  10%.  The  GPM  results  consistently  lower  then  those  of  the  Gaussian 
results  by  approximately  5%  and  this  is  true  for  all  of  the  trajectories  analyzed. 

The  difference  between  the  GAMESS  frequencies  and  those  calculated  by  GPM  can 
be  attributed  to  two  contributing  effects.  The  first  can  be  explained  by  observing  that 
the  frequencies  obtained  from  SCF  calculations  are  single  point  calculations  assuming 
a quadratic  potential  for  the  equilibrium.  The  GPM  analyzes  actual  nuclear  motion 
which  experience  anharmonic  effects  due  to  a breakdown  of  the  quadratic  approxima- 
tion as  bonds  are  stretched  from  equilibrium.  These  anharmonic  effects  tend  to  lower 
the  fundamental  frequency  to  first  order.  The  second  contributing  effect  comes  from 
an  understanding  of  the  implementation  of  the  END  in  the  ENDyne  code.  In  this  cur- 
rent implementation  the  nuclei  move  dynamically  through  an  average  potential  energy 
surface  which  for  the  process  under  investigation  was  a superposition  of  the  H20  and 
H20+  states.  The  probability  for  charge  transfer  determines  the  proportion  of  this  super- 
position. Under  the  influence  of  this  average  state  the  bond  strength  between  the  nuclei 
would  be  decreased  leading  to  a lower  frequency  for  the  normal  modes. 

The  unique  ability  of  ENDyne  to  calculate  the  full  dynamics  of  molecular  collisions 
allows  observation  of  molecular  processes  not  previously  possible  in  either  theory  or 
experiment.  The  experimental  results  are  far  more  statistical  in  nature  and  in  order  to 
compare  with  those  results  many  ENDyne  trajectories  must  be  analyzed. 

The  state  resolved  DCS  were  calculated  using  GPM  and  the  quasiclassical  coherent 
states.  As  we  observe  that  in  Fig.  7.17  and  Fig.  7.18  the  qualitative  structure  of  the  theo- 
retical and  experimental  state  resolved  DCS  are  in  strong  agreement.  The  comparison  of 
the  Non-Transfer  (NT)  DCS  shown  Fig.  7.19  exhibits  excellent  agreement  over  a broad 
range  of  deflection  angles.  We  note  here  that  all  the  experimental  results  used  arbitrary 
units,  we  normalized  the  results  to  5°  the  location  of  the  theoretical  and  experimental 
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V\ 

V 2 

K3 

freq  ( cm  x) 

3562 

1595 

3756 

STO  - 3G 

4140 

2170 

4390 

Prony 

4014 

2093 

4299 

E vib  (3-U) 

.009 

.0015 

.0024 

zy2  — Euik. 
U ~ hu; 

0.44 

0.16 

0.13 

Table  7.6: 


rainbow  angle.  The  ground  state  DCS  in  Fig  7.20  also  shows  excellent  agreement  for 
lower  deflection  angles  starting  to  diverge  at  higher  angles.  The  agreement  continues  to 
be  excellent  for  excitation  of  bending  mode  in  Fig.  7.21.  The  DCS  for  the  stretching 
modes  are  not  quite  as  good  as  seen  in  Fig.  7.22.  We  see  in  Fig.  7.23  that  the  asymmetric 
stretching  mode  makes  little  or  no  contribution  to  the  overall  stretching  DCS. 

In  Fig.  7.23  answers  the  inquiry  made  by  Toennies,  et  al.  by  resolving  the  asym- 
metric and  symmetric  stretching  modes.  In  their  paper  they  argue  that  the  during  the 
collision  the  presence  of  the  H + H20+  quasimolecular  state  would  lead  to  a preferred 
excitation  of  the  symmetric  stretching  mode.  The  resolution  of  these  modes  was  not 
possible  in  their  experiment  due  to  the  small  energy  gap  between  the  symmetric  and 
asymmetric  stretching  modes.  Our  ability  to  follow  the  full  nonadiabatic  dynamics  of 
the  collision  in  our  theoretical  model  has  allowed  us  to  investigate  this  question  with 
certainty. 

The  most  interesting  feature  of  this  theoretical  analysis  can  be  observed  in  where 
the  asymmetric  bending  mode  has  a very  low  cross  section.  Toennies,  et  al. [57]  ex- 
periments were  unable  to  resolve  the  symmetric  and  asymmetric  bending  modes  due 
to  small  energy  gap  between  them.  They  hypothesized  the  presence  of  the  H -I-  H20+ 
quasimolecular  state  during  the  reaction  process  would  lead  to  an  enhancement  of  the 
symmetric  stretching  and  bending  modes.  As  we  have  observed  the  low  asymmetric 
SRDCS  points  strongly  to  the  presence  of  the  H + H20+  quasimolecular  state.  More- 
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Figure  7.17:  Vibrationally  resolved  DCS  for  expermintal  results:  Source  Toennies[57] 


over  these  results  can  be  attributed  to  the  ability  for  the  END  theory  and  the  ENDyne 
program  to  simulate  the  full  dynamics  of  a nonadiabatic[59]  process  (see  Fig.  7.25. 
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Figure  7.19:  Theoretical  non-transfer  differential  cross  section  for  the  non-transfer  col- 
lision H+  T-  H20  compared  with  experiment. 
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7.20:  Ground  State  resolved  DCS:  Comparison  of  Theory  and  Experiment 


Figure  7.21:  Bending  Mode  resolved  DCS:  Comparison  of  Theory  and  Experiment 
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Figure  7.22:  Combined  stretching  modes  resolved  DCS:  Comparison  of  Theory  and 
Experiment 
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Figure  7.23:  Asymmetric  and  symmetric  modes  resolved  DCS:  Comparison  of  Theory 
and  Experiment 
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Figure  7.24:  Comparison  of  the  geometric  nature  of  H20  and  H20+ 


J.  Phys.  Chem.  A,  Vol.  103,  No.  35,  1999  7117 


Figure  7.25:  Potential  Energy  Surfaces  for  the  collison  process  H+  + H20:  Source  Di 
Giacomo, et  al.[59] 


CHAPTER  8 

ROTATIONAL  CALCULATION  RESULTS 


As  we  have  stated,  long  after  the  collision  has  taken  place  and  the  participants  are  spa- 
tially separated  we  may  consider  the  angular  momentum  of  the  H20*  to  be  conserved. 
Unlike  the  vibrational  analysis  we  need  only  calculate  the  angular  momentum  once  per 
trajectory.  Using  the  classical  formula  for  the  angular  momentum  in  the  laboratory 
frame  we  have 


We  calculated  this  for  each  trajectory  of  each  orientation;  as  an  example  the  results  for 
one  typical  orientation  are  displayed  in  Fig.  8.1. 

The  distribution  of  the  angular  momentum  modulus  as  a function  of  impact  param- 
eter for  this  orientation  displays  some  interesting  structure.  It  mostly  rises  for  lower 
impact  parameters  until  it  drops  precipitously.  The  drop  can  be  attributed  to  the  impact- 
ing proton  striking  the  water  molecule  head-on  as  opposed  to  a grazing  blow  that  would 
result  in  imparting  angular  momentum  to  the  water  molecule.  The  modulus  of  the  an- 
gular momentum  drops  for  large  impact  parameters  but  does  not  go  to  zero  as  quickly  as 
the  deflection  function.  This  is  the  result  of  the  long  range  coulomb  interaction  between 
the  proton  and  the  strong  dipole  moment  of  the  water  molecule. 

Once  the  calculation  of  the  classical  angular  momentum  has  been  accomplished  we 
may  analyze  its  physical  parameters  and  calculate  the  probability  to  be  in  a given  state 
associated  with  the  angular  momentum  operator  by  using 


(8.1) 
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(8.2) 


[(21)  !]2 


pd+M)(i  _p)(I-M)  (8  3) 


x</(I  + K)(l  — q)^  K)  exp (— r) 


I!  (I  + M)!(I-M)!(I  + K)!(I-K)! 
q( I + K)(l  - q)(1~K)  exp(-r)^, 


(8.4) 


from  Eqn.  5.47.  The  physical  parameters  we  need  are  /3  the  polar  angle,  9 the  azimuthal 
angle,  and  r the  modulus  of  the  body  frame  classical  angular  momentum.  These  param- 
eters are  given  by 


We  restricted  ourselves  to  an  analysis  of  the  angular  momentum  quantum  numbers.  This 
gives  a Poisson’s  distribution  centered  on  the  angular  momentum  modulus  given  by  the 
classical  value  for  the  magnitude  of  the  angular  momentum.  As  an  example  we  show 
the  distribution  associated  with  a particular  trajectory  in  Fig.  8.2  which  displays  the 
expected  Poisson’s  distribution. 

Using  the  same  body  of  data  used  in  our  vibrational  analysis  we  calculated  an  ex- 
perimentally relevant  parameter  the  integral  cross  section  given  by 


r 


x 


(8.5) 


(8.6) 
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We  employed  the  probability  factor  given  by  Eqn.  5.47  to  calculate  the  rotationally 
resolved  integral  cross  section  for  each  orientation  of  the  water  target.  Employing  the 
scheme  for  rotational  averaging  given  by  Jacquemin,  et  al.  [56]  we  were  able  to  generate 
Fig.  8.3.  We  see  that  resolved  states  show  very  little  in  the  way  of  interesting  features. 
The  large  hump  centered  at  approximately  \L\  — 3 au  is  very  similar  to  a Poisson’s  dis- 
tribution. The  rotational  analysis  reveals  that  the  rotational  excitation  for  this  collision 
at  this  energy  is  rather  low. 

In  light  of  the  analysis  in  Ch.  5 the  rotational  CS  are  only  approximately  accurate 
and  only  quasi-classical  in  the  high  angular  momentum  limit.  This  suggests  that  this 
set  of  data  is  not  a good  candidate  for  our  rotational  analysis  method.  We  regard  the 
preceding  analysis  and  an  illustrative  exercise.  Future  work  is  intended  to  explore  and 
develop  this  method  for  rotations  and  compare  to  several  experimental  results.  Possi- 
ble experimental  work  we  intend  to  investigate  include  the  charge  transfer  process  of 
Ar+  + O2  — > Ar  + by  Ng,  et  al.  [60]. 
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Deflection  Function 


Angular  Momentum  Versus  Impact  Parameter 
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Figure  8.1:  Comparison  of  the  deflection  function  and  classical  angular  momentum  as  a 
function  of  impact  parameter 
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Distribuition  for  Ang.  Quantum  Number 
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Figure  8.2:  Distribution  of  rotational  coherent  states  over  angular  momentum  eigenval- 
ues 


Rotationally  Resolved  Integral  Cross  Section 
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Figure  8.3:  Rotationally  resolved  integral  cross  section  as  a function  of  angular  momen- 
tum eigenvalues 


CHAPTER  9 
CONCLUSION 


Quantum  chemical  ab  initio  calculations  are  often  quite  computationally  expensive  and 
are  often  simplified  by  keeping  the  nuclei  classical.  This  approximation  is  valid  for 
a broad  range  of  energies  and  chemical  processes.  Exploiting  the  connection  between 
quasiclassical  coherent  states  and  semiclassical  dynamics  we  were  able  to  derive  numer- 
ical techniques  to  map  quantum  mechanical  states  to  classical  trajectories.  We  restate 
the  overriding  principle  of  this  dissertation:  Using  coherent  states  and  a posteriori  analy- 
sis we  can  extract  state-resolved  results  from  a semiclassical  description  using  classical 
nuclei.  Specifically,  we  strove  to  study  ro vibrational  dynamics  and  the  state-resolved 
differential  cross  sections  in  molecular  scattering.  Building  a technique  on  three  sepa- 
rate pillars — Coherent  States,  analysis  of  classical  dynamics,  and  molecular  trajectories 
calculated  using  QCSD  END  theory — we  successfully  studied  state  resolved  dynamics 
of  several  interesting  systems. 

We  applied  our  technique  to  the  theoretical  examination  of  crossed  beam  scattering 
experiments.  Using  the  Electron  Nuclear  Dynamics  theory  to  calculate  the  dynamics  we 
simulated  experimental  conditions.  In  the  theory’s  first  level  of  approximation  the  nuclei 
are  treated  classically  which  is  implemented  in  the  ENDyne  program.  The  Generalized 
Prony’s  method  we  developed  and  described  was  the  necessary  tool  to  successfully  link 
the  quasiclassical  coherent  states  with  the  classical  description  of  the  nuclei  and  compare 
to  state  resolved  experimental  data.  In  our  analysis  of  vibrationally  resolved  differential 
cross  sections  of  inelastic  scattering  of  H+  + H20  at  46  eV  we  were  able  to  successfully 
relate  vibrational  coherent  states  with  the  motion  of  the  classical  nuclei.  Moreover, 
we  were  able  to  investigate  the  nonadiabatic  mechanisms  involved  in  this  scattering 
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experiment  and  add  theoretical  insight  into  the  presence  of  intermediate  quasi  molecular 
states. 

Our  thorough  investigation  into  the  effectiveness  of  GPM  in  analyzing  nuclear  vibra- 
tional trajectories  and  favorable  results  obtained  from  utilizing  these  results  allows  us  to 
conclude  that  our  theoretical  technique  to  describe  classical  nuclei  with  quasiclassical 
coherent  states  is  a powerful  tool.  This  technique  expands  the  ability  of  the  researcher 
to  analyze  and  interpret  their  data  and  compare  to  experiment.  We  expect  these  devel- 
opments to  find  broad  application  in  state  resolved  reaction  dynamics  since  it  creates 
the  opportunity  to  gain  insight  into  theoretical  calculations  without  using  full  quantum 
nuclei  significantly  reducing  computation  time. 
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